Square form factorization

Daniel Shanks's Square Form Factorization (SquFoF).

Task
Square form factorization
You are encouraged to solve this task according to the task description, using any language you may know.
Task.

Invented around 1975, ‘On a 32-bit computer, SquFoF is the clear champion factoring algorithm for numbers between 1010 and 1018, and will likely remain so.’


An integral binary quadratic form is a polynomial f(x,y) = ax2 + bxy + cy2 with integer coefficients and discriminant D = b24ac. For each positive discriminant there are multiple forms (a, b, c).

The next form in a periodic sequence (cycle) of adjacent forms is found by applying a reduction operator rho, essentially a variant of Euclid's algorithm for finding the continued fraction of a square root. Using floor(√N), rho constructs a principal form (1, b, c) with D = 4N.

SquFoF is based on the existence of cycles containing ambiguous forms, with the property that a divides b. They come in pairs of associated forms (a, b, c) and (c, b, a) called symmetry points. If an ambiguous form is found (there is one for each divisor of D), write the discriminant as (ak)24ac = a(a·k24c) = 4N and (if a is not equal to 1 or 2) N is split.

Shanks used square forms to jump to a random ambiguous cycle. Fact: if any form in an ambiguous cycle is squared, that square form will always land in the principal cycle. Conversely, the square root of any form in the principal cycle lies in an ambiguous cycle. (Possibly the principal cycle itself).

A square form is easy to find: the last coefficient c is a perfect square. This happens about once every ∜N-th cycle step and for even indices only. Let rho compute the inverse square root form and track the ambiguous cycle backward until the symmetry point is reached. (Taking the inverse reverses the cycle). Then a or a/2 divides D and therefore N.

To avoid trivial factorizations, Shanks created a list (queue) to hold small coefficients appearing early in the principal cycle, that may be roots of square forms found later on. If these forms are skipped, no roots land in the principal cycle itself and cases a = 1 or a = 2 do not happen.

Sometimes the cycle length is too short to find a proper square form. This is fixed by running five instances of SquFoF in parallel, with input N and 3, 5, 7, 11 times N; the discriminants then will have different periods. If N is prime or the cube of a prime, there are improper squares only and the program will duly report failure.

Reference.

[1] A detailed analysis of SquFoF (2007)


ALGOL 68

Works with: ALGOL 68G version Any - tested with release 2.8.3.win32

Assumes LONG INT is long enough - this is true in ALGOL 68G versioins 2 and 3.

Translation of: Wren – ...but only using a single size of integer
BEGIN # Daniel Shanks's Square Form Factorization (SquFoF) - based on the Wren sample #

    MODE INTEGER = LONG INT;                                  # large enough INT type #
    PROC(LONG REAL)LONG REAL size sqrt = long sqrt;         # sqrt for INTEGER values #

    []INTEGER multipliers = ( 1, 3, 5, 7, 11, 3 * 5, 3 * 7, 3 * 11
                            , 5 * 7, 5 * 11, 7 * 11, 3 * 5 * 7, 3 * 5 * 11
                            , 3 * 7 * 11, 5 * 7 * 11, 3 * 5 * 7 * 11
                            );
    PROC gcd = ( INTEGER x, y )INTEGER:                               # iterative gcd #
         BEGIN
            INTEGER a := x, b := y;
            WHILE b /= 0 DO
               INTEGER next a = b;
               b := a MOD b;
               a := next a
            OD;
            ABS a
         END # gcd # ;

    PROC squfof = ( INTEGER n )INTEGER:
         IF  INTEGER s = ENTIER ( size sqrt( n ) + 0.5 );
             s * s = n
         THEN s
         ELSE INTEGER result := 0;
              FOR multiplier FROM LWB multipliers TO UPB multipliers WHILE result = 0 DO
                  INTEGER d       = n * multipliers[ multiplier ];
                  INTEGER pp     := ENTIER size sqrt( d );
                  INTEGER p prev := pp;
                  INTEGER po      = p prev;
                  INTEGER q prev := 1;
                  INTEGER qq     := d - ( po * po );
                  INTEGER l       = ENTIER size sqrt( s * 8 );
                  INTEGER bb      = 3 * l;
                  INTEGER i      := 2;
                  INTEGER b      := 0;
                  INTEGER q      := 0;
                  INTEGER r      := 0;
                  BOOL    again  := TRUE;
                  WHILE i < bb AND again DO
                      b  := ( po + pp ) OVER qq;
                      pp := ( b * qq ) - pp;
                      q  := qq;
                      qq := q prev + ( b * ( p prev - pp ) );
                      r  := ENTIER ( size sqrt( qq ) + 0.5 );
                      IF i MOD 2 = 0 THEN again := r * r /= qq FI;
                      IF again THEN
                          q prev := q;
                          p prev := pp;
                          i     +:= 1
                      FI
                  OD;
                  IF i < bb THEN
                      b      := ( po - pp ) OVER r;
                      p prev := pp := ( b * r ) + pp;
                      q prev := r;
                      qq     := ( d - ( p prev * p prev ) ) OVER q prev;
                      i      := 0;
                      WHILE
                          b      := ( po + pp ) OVER qq;
                          p prev := pp;
                          pp     := ( b * qq ) - pp;
                          q      := qq;
                          qq     := q prev + ( b * ( p prev - pp ) );
                          q prev := q;
                          i     +:= 1;
                          pp /= p prev
                      DO SKIP OD
                  FI;
                  r := gcd( n, q prev );
                  IF r /= 1 AND r /=n THEN result := r FI
              OD;
              result
         FI # squfof # ;

    []INTEGER examples = (                2501,               12851
                         ,               13289,               75301
                         ,              120787,              967009
                         ,              997417,             7091569
                         ,            13290059,            42854447
                         ,           223553581,          2027651281
                         ,         11111111111,        100895598169
                         ,       1002742628021,      60012462237239
                         ,     287129523414791,    9007199254740931
                         ,   11111111111111111,  314159265358979323
                         ,  384307168202281507,  419244183493398773
                         ,  658812288346769681,  922337203685477563
                         , 1000000000000000127, 1152921505680588799
                         , 1537228672809128917, 4611686018427387877
                     );

    print( ( "Integer                  Factor Quotient",   newline ) );
    print( ( "----------------------------------------", newline ) );
    FOR example FROM LWB examples TO UPB examples DO
        INTEGER n    = examples[ example ];
        INTEGER fact = squfof( n );
        STRING  quot = IF fact = 0 THEN "fail" ELSE whole( n OVER fact, 0 ) FI;
        print( ( whole( n, -20 ), " ", whole( fact, -10 ), " ", quot, newline ) )
    OD
END
Output:
Integer                  Factor Quotient
----------------------------------------
                2501         61 41
               12851         71 181
               13289        137 97
               75301        293 257
              120787         43 2809
              967009       1609 601
              997417        257 3881
             7091569       2663 2663
            13290059       3119 4261
            42854447       9689 4423
           223553581      11213 19937
          2027651281      46061 44021
         11111111111      21649 513239
        100895598169     112303 898423
       1002742628021          0 fail
      60012462237239    6862753 8744663
     287129523414791    6059887 47381993
    9007199254740931   10624181 847801751
   11111111111111111    2071723 5363222357
  314159265358979323  317213509 990371647
  384307168202281507  415718707 924440401
  419244183493398773   48009977 8732438749
  658812288346769681   62222119 10588072199
  922337203685477563  110075821 8379108103
 1000000000000000127  111756107 8948056861
 1152921505680588799  139001459 8294312261
 1537228672809128917   26675843 57626245319
 4611686018427387877  343242169 13435662733

C

Based on Wp

From the Wikipedia entry for Shanks's square forms factorization [[2]]

#include <math.h>
#include <stdio.h>

#define nelems(x) (sizeof(x) / sizeof((x)[0]))

const unsigned long multiplier[] = {1, 3, 5, 7, 11, 3*5, 3*7, 3*11, 5*7, 5*11, 7*11, 3*5*7, 3*5*11, 3*7*11, 5*7*11, 3*5*7*11};

unsigned long long gcd(unsigned long long a, unsigned long long b)
{
    while (b != 0)
    {
        a %= b;
        a ^= b;
        b ^= a;
        a ^= b;
    }

    return a;
}

unsigned long long SQUFOF( unsigned long long N )
{
    unsigned long long D, Po, P, Pprev, Q, Qprev, q, b, r, s;
    unsigned long L, B, i;
    s = (unsigned long long)(sqrtl(N)+0.5);
    if (s*s == N) return s;
    for (int k = 0; k < nelems(multiplier) && N <= 0xffffffffffffffff/multiplier[k]; k++) {
        D = multiplier[k]*N;
        Po = Pprev = P = sqrtl(D);
        Qprev = 1;
        Q = D - Po*Po;
        L = 2 * sqrtl( 2*s );
        B = 3 * L;
        for (i = 2 ; i < B ; i++) {
            b = (unsigned long long)((Po + P)/Q);
            P = b*Q - P;
            q = Q;
            Q = Qprev + b*(Pprev - P);
            r = (unsigned long long)(sqrtl(Q)+0.5);
            if (!(i & 1) && r*r == Q) break;
            Qprev = q;
            Pprev = P;
        };
        if (i >= B) continue;
        b = (unsigned long long)((Po - P)/r);
        Pprev = P = b*r + P;
        Qprev = r;
        Q = (D - Pprev*Pprev)/Qprev;
        i = 0;
        do {
            b = (unsigned long long)((Po + P)/Q);
            Pprev = P;
            P = b*Q - P;
            q = Q;
            Q = Qprev + b*(Pprev - P);
            Qprev = q;
            i++;
        } while (P != Pprev);
        r = gcd(N, Qprev);
        if (r != 1 && r != N) return r;
    }
    return 0;
}

int main(int argc, char *argv[]) {
    int i;
    const unsigned long long data[] = {
        2501,
        12851,
        13289,
        75301,
        120787,
        967009,
        997417,
        7091569,
        13290059,
        42854447,
        223553581,
        2027651281,
        11111111111,
        100895598169,
        1002742628021,
        60012462237239,
        287129523414791,
        9007199254740931,
        11111111111111111,
        314159265358979323,
        384307168202281507,
        419244183493398773,
        658812288346769681,
        922337203685477563,
        1000000000000000127,
        1152921505680588799,
        1537228672809128917,
        4611686018427387877};

    for(int i = 0; i < nelems(data); i++) {
        unsigned long long example, factor, quotient;
        example = data[i];
        factor = SQUFOF(example);
        if(factor == 0) {
            printf("%llu was not factored.\n", example);
        }
        else {
            quotient = example / factor;
            printf("Integer %llu has factors %llu and %llu\n",
               example, factor, quotient); 
        }
    }
}
Output:
Integer 2501 has factors 61 and 41
Integer 12851 has factors 71 and 181
Integer 13289 has factors 137 and 97
Integer 75301 has factors 293 and 257
Integer 120787 has factors 43 and 2809
Integer 967009 has factors 1609 and 601
Integer 997417 has factors 257 and 3881
Integer 7091569 has factors 2663 and 2663
Integer 13290059 has factors 3119 and 4261
Integer 42854447 has factors 9689 and 4423
Integer 223553581 has factors 11213 and 19937
Integer 2027651281 has factors 46061 and 44021
Integer 11111111111 has factors 21649 and 513239
Integer 100895598169 has factors 112303 and 898423
1002742628021 was not factored.
Integer 60012462237239 has factors 6862753 and 8744663
Integer 287129523414791 has factors 6059887 and 47381993
Integer 9007199254740931 has factors 10624181 and 847801751
Integer 11111111111111111 has factors 2071723 and 5363222357
Integer 314159265358979323 has factors 317213509 and 990371647
Integer 384307168202281507 has factors 415718707 and 924440401
Integer 419244183493398773 has factors 48009977 and 8732438749
Integer 658812288346769681 has factors 62222119 and 10588072199
Integer 922337203685477563 has factors 110075821 and 8379108103
1000000000000000127 was not factored.
Integer 1152921505680588799 has factors 139001459 and 8294312261
1537228672809128917 was not factored.
Integer 4611686018427387877 has factors 343242169 and 13435662733

Classical heuristic

See Discussion.

//SquFoF: minimalistic version without queue.
//Classical heuristic. Tested: tcc 0.9.27
#include <math.h>
#include <stdio.h>

//input maximum
#define MxN ((unsigned long long) 1 << 62)

//reduce indefinite form
#define rho(a, b, c) {        \
  t = c; c = a; a = t; t = b; \
   q = (rN + b) / a;          \
   b = q * a - b;             \
   c += q * (t - b); }

//initialize
#define rhoin(a, b, c) {      \
   rho(a, b, c)  h = b;       \
   c = (mN - h * h) / a; }

#define gcd(a, b) while (b) { \
   t = a % b; a = b; b = t; }

//multipliers
const unsigned long m[] = {1, 3, 5, 7, 11, 0};

//square form factorization
unsigned long squfof( unsigned long long N ) {
unsigned long a, b, c, u, v, w, rN, q, t, r;
unsigned long long mN, h;
int i, ix, k = 0;

   if ((N & 1)==0) return 2;

   h = floor(sqrt(N)+ 0.5);
   if (h * h == N) return h;

   while (m[k]) {
      if (k && N % m[k]==0) return m[k];
      //check overflow m * N
      if (N > MxN / m[k]) break;
      mN = N * m[k++];

      r = floor(sqrt(mN));
      h = r; //float64 fix
      if (h * h > mN) r -= 1;
      rN = r;

      //principal form
      b = r; c = 1;
      rhoin(a, b, c)

      //iteration bound
      ix = floor(sqrt(2*r)) * 4;

      //search principal cycle
      for (i = 2; i < ix; i += 2) {
         rho(a, b, c)
         //even step

         r = floor(sqrt(c)+ 0.5);
         if (r * r == c) {
            //square form found

            //inverse square root
            v = -b; w = r;
            rhoin(u, v, w)

            //search ambiguous cycle
            do { r = v;
              rho(u, v, w)
            } while (v != r);
            //symmetry point

            h = N; gcd(h, u)
            if (h != 1) return h;
         }
         rho(a, b, c)
         //odd step
      }
   }
   return 1;
}

void main(void) {
const unsigned long long data[] = {
   2501,
   12851,
   13289,
   75301,
   120787,
   967009,
   997417,
   7091569,

   5214317,
   20834839,
   23515517,
   33409583,
   44524219,

   13290059,
   223553581,
   2027651281,
   11111111111,
   100895598169,
   1002742628021,
   60012462237239,
   287129523414791,
   9007199254740931,
   11111111111111111,
   314159265358979323,
   384307168202281507,
   419244183493398773,
   658812288346769681,
   922337203685477563,
   1000000000000000127,
   1152921505680588799,
   1537228672809128917,
   4611686018427387877,
   0};

   unsigned long long N, f;
   int i = 0;

   while (1) {
      N = data[i++];
      //scanf("%llu", &N);
      if (N < 2) break;

      printf("N = %llu\n", N);

      f = squfof(N);
      if (N % f) f = 1;

      if (f == 1) printf("fail\n\n");
      else printf("f = %llu  N/f = %llu\n\n", f, N/f);
   }
}
Showing problem cases only:
...
N = 5214317
f = 73  N/f = 71429

N = 20834839
f = 3361  N/f = 6199

N = 23515517
f = 53  N/f = 443689

N = 33409583
f = 991  N/f = 33713

N = 44524219
f = 593  N/f = 75083
...
N = 1000000000000000127
f = 111756107  N/f = 8948056861
...
N = 1537228672809128917
f = 26675843  N/f = 57626245319
...

C++

#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <random>

uint64_t test_value = 0;
uint64_t sqrt_test_value = 0;

class BQF { // Binary quadratic form
public:
	BQF(const uint64_t& a, const uint64_t& b, const uint64_t& c) : a(a), b(b), c(c) {
		q = ( sqrt_test_value + b ) / c;
		bb = q * c - b;
	}

	BQF rho() {
		return BQF(c, bb, a  + q * ( b - bb ));
	}

	BQF rho_inverse() {
		return BQF(c, bb, ( test_value - bb * bb ) / c);
	}

	uint64_t a, b, c;
private:
	uint64_t q, bb;
};

uint64_t squfof(const uint64_t& number) {
	const uint32_t sqrt = std::sqrt(number);
	if ( sqrt * sqrt == number ) {
		return sqrt;
	}

	test_value = number;
	sqrt_test_value = std::sqrt(test_value);

	// Principal form
	BQF form(0, sqrt_test_value, 1);
	form = form.rho_inverse();

	// Search principal cycle
	for ( uint32_t i = 0; i < 4 * std::sqrt(2 * sqrt_test_value); i += 2 ) {
		// Even step
		form = form.rho();

		uint64_t sqrt_c = std::sqrt(form.c);
		if ( sqrt_c * sqrt_c == form.c ) { // Square form found
			// Inverse square root
			BQF form_inverse(0, -form.b, sqrt_c);
			form_inverse = form_inverse.rho_inverse();

			// Search ambiguous cycle
			uint64_t previous_b = 0;
			do {
				previous_b = form_inverse.b;
				form_inverse = form_inverse.rho();
			} while ( form_inverse.b != previous_b );

			// Symmetry point
			const uint64_t g = std::gcd(number, form_inverse.a);
			if ( g != 1 ) {
				return g;
			}
		}

		// Odd step
		form = form.rho();
	}

	if ( number % 2 == 0 ) {
		return 2;
	}
	return 0; // Failed to factorise, possibly a prime number
}

int main() {
	std::random_device random;
	std::mt19937 generator(random());
	const uint64_t lower_limit = 100'000'000'000'000'000;
	std::uniform_int_distribution<uint64_t> distribution(lower_limit, 10 * lower_limit);

	for ( uint32_t i = 0; i < 20; ++i ) {
		uint64_t test = distribution(random);
		uint64_t factor = squfof(test);

		if ( factor == 0 ) {
			std::cout << test << " - failed to factorise" << std::endl;
		} else {
			std::cout << test << " = " << factor << " * " << test / factor << std::endl;
		}
		std::cout << std::endl;
   }
}
Output:
822140815871714649 = 141 * 5830785928168189

473377979025428817 = 3 * 157792659675142939

482452941918160803 = 4410431 * 109389069213

165380937127655630 = 65438 * 2527292049385

191677853606692475 = 7589219 * 25256598025

480551815975206727 = 2843 * 169029833265989

178710207362206205 = 5 * 35742041472441241

484660189375949842 = 1094 * 443016626486243

758704390319635770 = 1605 * 472713015775474

820453356193182720 = 97280 * 8433936638499

706982627912630220 = 121273 * 5829678724140

614913973550671312 = 437204432 * 1406467841

601482456081568543 = 131 * 4591469130393653

610533314488947626 = 14 * 43609522463496259

336343281182924332 = 70108 * 4797502156429

308127213282933401 = 7 * 44018173326133343

582455924775519843 = 3 * 194151974925173281

694215100094443276 = 32070628 * 21646445467

398821795604697523 = 181 * 2203435334832583

477964959783291032 = 517029608 * 924444079

EasyLang

Translation of: C
multiplier[] = [ 1 3 5 7 11 3 * 5 3 * 7 3 * 11 5 * 7 5 * 11 7 * 11 3 * 5 * 7 3 * 5 * 11 3 * 7 * 11 5 * 7 * 11 3 * 5 * 7 * 11 ]
func gcd a b .
   while b <> 0
      a = a mod b
      swap a b
   .
   return a
.
func squfof N .
   s = floor (sqrt N + 0.5)
   if s * s = N
      return s
   .
   for multiplier in multiplier[]
      if N > 9007199254740992 / multiplier
         print "Number " & N & " is too big"
         break 1
      .
      D = multiplier * N
      P = floor sqrt D
      Po = P
      Pprev = P
      Qprev = 1
      Q = D - Po * Po
      L = 2 * floor sqrt (2 * s)
      B = 3 * L
      for i = 2 to B - 1
         b = (Po + P) div Q
         P = b * Q - P
         q = Q
         Q = Qprev + b * (Pprev - P)
         r = floor (sqrt Q + 0.5)
         if i mod 2 = 0 and r * r = Q
            break 1
         .
         Qprev = q
         Pprev = P
      .
      if i < B
         b = (Po - P) div r
         P = b * r + P
         Pprev = P
         Qprev = r
         Q = (D - Pprev * Pprev) / Qprev
         i = 0
         repeat
            b = (Po + P) div Q
            Pprev = P
            P = b * Q - P
            q = Q
            Q = Qprev + b * (Pprev - P)
            Qprev = q
            i += 1
            until P = Pprev
         .
         r = gcd N Qprev
         if r <> 1 and r <> N
            return r
         .
      .
   .
   return 0
.
data[] = [ 2501 12851 13289 75301 120787 967009 997417 7091569 13290059 42854447 223553581 2027651281 11111111111 100895598169 1002742628021 60012462237239 287129523414791 9007199254740931 ]
for example in data[]
   factor = squfof example
   if factor = 0
      print example & " was not factored."
   else
      quotient = example / factor
      print example & " has factors " & factor & " " & quotient
   .
.

FreeBASIC

' ***********************************************
'subject: Shanks's square form factorization:
'         ambiguous forms of discriminant 4N
'         give factors of N.
'tested : FreeBasic 1.08.1


'------------------------------------------------
const MxN = culngint(1) shl 62
'input maximum

const qx = (1 shl 5) - 1
'queue size

type arg
   'squfof arguments
   as ulong m, f
   as integer vb
end type

type bqf
   declare sub rho ()
   'reduce indefinite form
   declare function issq (byref r as ulong) as integer
   'return -1 if c is square, set r:= sqrt(c)
   declare sub qform (byref g as string, byval t as integer)
   'print binary quadratic form #t (a, 2b, c)

   as ulong rN, a, b, c
   as integer vb
end type

type queue
   declare sub enq (byref P as bqf)
   'enqueue P.c, P.b if appropriate
   declare function pro (byref P as bqf, byval r as ulong) as integer
   'return -1 if a proper square form is found

   as ulong a(qx), L, m
   as integer k, t
end type

'global variables
dim shared N as ulongint
'the number to split

dim shared flag as integer
'signal to end all threads

dim shared as ubyte q1024(1023), q3465(3464)
'quadratic residue tables


'------------------------------------------------
sub bqf.rho ()
dim as ulong q, t
   swap a, c
   'residue
   q = culng(rN + b) \ a
   t = b: b = q * a - b
   'pseudo-square
   c += q * (t - b)
end sub

'initialize form
#macro rhoin(F)
   F.rho : h = F.b
   F.c = (mN - h * h) \ F.a
#endmacro

function bqf.issq (byref r as ulong) as integer
if q1024(c and 1023) andalso q3465(c mod 3465) then
   '98.6% non-squares filtered
   r = culng(sqr(c))
   if r * r = c then return -1
end if
issq = 0
end function

sub bqf.qform (byref g as string, byval t as integer)
if vb = 0 then exit sub
dim as longint u = a, v = b, w = c
   if t and 1 then
      w = -w
   else
      u = -u
   end if
   v shl= 1
   print g;str(t);" = (";u;",";v;",";w;")"
end sub

'------------------------------------------------
#macro red(r, a)
   r = iif(a and 1, a, a shr 1)
   if m > 2 then
      r = iif(r mod m, r, r \ m)
   end if
#endmacro

sub queue.enq (byref P as bqf)
dim s as ulong
   red(s, P.c)
   if s < L then
      'circular queue
      k = (k + 2) and qx
      if k > t then t = k
      'enqueue P.b, P.c
      a(k) = P.b mod s
      a(k + 1) = s
   end if
end sub

function queue.pro (byref P as bqf, byval r as ulong) as integer
dim as integer i, sw
   'skip improper square forms
   for i = 0 to t step 2
      sw = (P.b - a(i)) mod r = 0
      sw and= a(i + 1) = r
      if sw then return 0
   next i
pro = -1
end function

'------------------------------------------------
sub squfof (byval ap as any ptr)
dim as arg ptr rp = cptr(arg ptr, ap)
dim as ulong L2, m, r, t, f = 1
dim as integer ix, i, j
dim as ulongint mN, h
'principal and ambiguous cycles
dim as bqf P, A
dim Q as queue

if (N and 1) = 0 then
   rp->f = 2 ' even N
   flag =-1: exit sub
end if

h = culngint(sqr(N))
if h * h = N then
   'N is square
   rp->f = culng(h)
   flag =-1: exit sub
end if

rp->f = 1
'multiplier
m = rp->m
if m > 1 then
   if (N mod m) = 0 then
      rp->f = m  ' m | N
      flag =-1: exit sub
   end if

   'check overflow m * N
   if N > (MxN \ m) then exit sub
end if
mN = N * m

r = int(sqr(mN))
'float64 fix
if culngint(r) * r > mN then r -= 1
P.rN = r
A.rN = r

P.vb = rp->vb
A.vb = rp->vb
'verbosity switch
if P.vb then print "r = "; r

Q.k = -2: Q.t = -1: Q.m = m
'Queue entry bounds
Q.L = int(sqr(r * 2))
L2 = Q.L * m shl 1

'principal form
P.b = r: P.c = 1
rhoin(P)
P.qform("P", 1)

ix = Q.L shl 2
for i = 2 to ix
   'search principal cycle

   if P.c < L2 then Q.enq(P)

   P.rho
   if (i and 1) = 0 andalso P.issq(r) then
      'square form found

      if Q.pro(P, r) then

         P.qform("P", i)
         'inverse square root
         A.b =-P.b: A.c = r
         rhoin(A): j = 1
         A.qform("A", j)

         do
            'search ambiguous cycle
            t = A.b
            A.rho: j += 1

            if A.b = t then
               'symmetry point
               A.qform("A", j)
               red(f, A.a)
               if f = 1 then exit do

               flag = -1
               'factor found
            end if
         loop until flag

      end if ' proper square
   end if ' square form

   if flag then exit for
next i

rp->f = f
end sub

'------------------------------------------------
data 2501
data 12851
data 13289
data 75301
data 120787
data 967009
data 997417
data 7091569
data 13290059
data 23515517
data 42854447
data 223553581
data 2027651281
data 11111111111
data 100895598169
data 1002742628021
data 60012462237239
data 287129523414791
data 9007199254740931
data 11111111111111111
data 314159265358979323
data 384307168202281507
data 419244183493398773
data 658812288346769681
data 922337203685477563
data 1000000000000000127
data 1152921505680588799
data 1537228672809128917
data 4611686018427387877
data 0

'main
'------------------------------------------------
const tx = 4
dim as double tim = timer
dim h(4) as any ptr
dim a(4) as arg
dim as ulongint f
dim as integer s, t

width 64, 30
cls

'tabulate quadratic residues
for t = 0 to 1540
   s = t * t
   q1024(s and 1023) =-1
   q3465(s mod 3465) =-1
next t

a(0).vb = 0
'set one verbosity switch only

a(0).m = 1
'multipliers
a(1).m = 3
a(2).m = 5
a(3).m = 7
a(4).m = 11

do
   print

   do : read N
   loop until N < MxN
   if N < 2 then exit do

   print "N = "; N

   flag = 0

   for t = 1 to tx + 1 step 2
      if t < tx then
         h(t) = threadcreate(@squfof, @a(t))
      end if

      squfof(@a(t - 1))
      f = a(t - 1).f

      if t < tx then
         threadwait(h(t))
         if f = 1 then f = a(t).f
      end if

      if f > 1 then exit for
   next t

   'assert
   if N mod f then f = 1

   if f = 1 then
      print "fail"
   else
      print "f = ";f;"  N/f = ";N \ f
   end if
loop

print "total time:"; csng(timer - tim); " s"
end
Examples:
N = 2501
f = 61  N/f = 41

N = 12851
f = 71  N/f = 181

N = 13289
f = 97  N/f = 137

N = 75301
f = 293  N/f = 257

N = 120787
f = 43  N/f = 2809

N = 967009
f = 1609  N/f = 601

N = 997417
f = 257  N/f = 3881

N = 7091569
f = 2663  N/f = 2663

N = 13290059
f = 3119  N/f = 4261

N = 23515517
f = 53  N/f = 443689

N = 42854447
f = 4423  N/f = 9689

N = 223553581
f = 11213  N/f = 19937

N = 2027651281
f = 44021  N/f = 46061

N = 11111111111
f = 21649  N/f = 513239

N = 100895598169
f = 112303  N/f = 898423

N = 1002742628021
fail

N = 60012462237239
f = 6862753  N/f = 8744663

N = 287129523414791
f = 6059887  N/f = 47381993

N = 9007199254740931
f = 10624181  N/f = 847801751

N = 11111111111111111
f = 2071723  N/f = 5363222357

N = 314159265358979323
f = 317213509  N/f = 990371647

N = 384307168202281507
f = 415718707  N/f = 924440401

N = 419244183493398773
f = 48009977  N/f = 8732438749

N = 658812288346769681
f = 62222119  N/f = 10588072199

N = 922337203685477563
f = 110075821  N/f = 8379108103

N = 1000000000000000127
f = 111756107  N/f = 8948056861

N = 1152921505680588799
f = 139001459  N/f = 8294312261

N = 1537228672809128917
f = 26675843  N/f = 57626245319

N = 4611686018427387877
f = 343242169  N/f = 13435662733

total time: 0.0170462 s

Go

Yet another solution based on the Wikipedia C code.

Rather than juggling with big.Int, I've just allowed the two 'awkward' cases to fail.

package main

import (
    "fmt"
    "math"
)

func isqrt(x uint64) uint64 {
    x0 := x >> 1
    x1 := (x0 + x/x0) >> 1
    for x1 < x0 {
        x0 = x1
        x1 = (x0 + x/x0) >> 1
    }
    return x0
}

func gcd(x, y uint64) uint64 {
    for y != 0 {
        x, y = y, x%y
    }
    return x
}

var multiplier = []uint64{
    1, 3, 5, 7, 11, 3 * 5, 3 * 7, 3 * 11, 5 * 7, 5 * 11, 7 * 11, 3 * 5 * 7, 3 * 5 * 11, 3 * 7 * 11, 5 * 7 * 11, 3 * 5 * 7 * 11,
}

func squfof(N uint64) uint64 {
    s := uint64(math.Sqrt(float64(N)) + 0.5)
    if s*s == N {
        return s
    }
    for k := 0; k < len(multiplier) && N <= math.MaxUint64/multiplier[k]; k++ {
        D := multiplier[k] * N
        P := isqrt(D)
        Pprev := P
        Po := Pprev
        Qprev := uint64(1)
        Q := D - Po*Po
        L := uint32(isqrt(8 * s))
        B := 3 * L
        i := uint32(2)
        var b, q, r uint64
        for ; i < B; i++ {
            b = uint64((Po + P) / Q)
            P = b*Q - P
            q = Q
            Q = Qprev + b*(Pprev-P)
            r = uint64(math.Sqrt(float64(Q)) + 0.5)
            if (i&1) == 0 && r*r == Q {
                break
            }
            Qprev = q
            Pprev = P
        }
        if i >= B {
            continue
        }
        b = uint64((Po - P) / r)
        P = b*r + P
        Pprev = P
        Qprev = r
        Q = (D - Pprev*Pprev) / Qprev
        i = 0
        for {
            b = uint64((Po + P) / Q)
            Pprev = P
            P = b*Q - P
            q = Q
            Q = Qprev + b*(Pprev-P)
            Qprev = q
            i++
            if P == Pprev {
                break
            }
        }
        r = gcd(N, Qprev)
        if r != 1 && r != N {
            return r
        }
    }
    return 0
}

func main() {
    examples := []uint64{
        2501,
        12851,
        13289,
        75301,
        120787,
        967009,
        997417,
        7091569,
        13290059,
        42854447,
        223553581,
        2027651281,
        11111111111,
        100895598169,
        1002742628021,
        60012462237239,
        287129523414791,
        9007199254740931,
        11111111111111111,
        314159265358979323,
        384307168202281507,
        419244183493398773,
        658812288346769681,
        922337203685477563,
        1000000000000000127,
        1152921505680588799,
        1537228672809128917,
        4611686018427387877,
    }
    fmt.Println("Integer              Factor     Quotient")
    fmt.Println("------------------------------------------")
    for _, N := range examples {
        fact := squfof(N)
        quot := "fail"
        if fact > 0 {
            quot = fmt.Sprintf("%d", N/fact)
        }
        fmt.Printf("%-20d %-10d %s\n", N, fact, quot)
    }
}
Output:
Integer              Factor     Quotient
------------------------------------------
2501                 61         41
12851                71         181
13289                137        97
75301                293        257
120787               43         2809
967009               1609       601
997417               257        3881
7091569              2663       2663
13290059             3119       4261
42854447             9689       4423
223553581            11213      19937
2027651281           46061      44021
11111111111          21649      513239
100895598169         112303     898423
1002742628021        0          fail
60012462237239       6862753    8744663
287129523414791      6059887    47381993
9007199254740931     10624181   847801751
11111111111111111    2071723    5363222357
314159265358979323   317213509  990371647
384307168202281507   415718707  924440401
419244183493398773   48009977   8732438749
658812288346769681   62222119   10588072199
922337203685477563   110075821  8379108103
1000000000000000127  0          fail
1152921505680588799  139001459  8294312261
1537228672809128917  0          fail
4611686018427387877  343242169  13435662733

J

Generally, this task should be accomplished in J using q:. Here we take an approach that's more comparable with the other examples on this page.

J does not have an unsigned fixed width integer type, which is one of the reasons that (in J) this algorithm is less optimal than advertised:

sqff=: {{
  s=. <.%:y 
  if. y=*:s do. s return. end.
  for_D. (x:y)*/:~*/@>,{1,each}.p:i.5 do.
    if. -.'integer'-:datatype D=. x:inv D do. break. end.
    P=. <.%:D
    Q=. 1, D-P*P
    lim=. <:6*<.%:2*s
    for_i. }.i.lim do.
      b=. <.(+/0 _1{P)%{:Q
      P=. P,|(b*{:Q)-{:P
      Q=. Q,|(_2{Q)+b*-/_2{.P
      if. 2|i do. if. (=<.&.%:){:Q do. break. end. end.
    end.
    if. i>:lim do. continue. end.
    Q=. <.%:{:Q
    b=. <.(-/0 _1{P)%Q
    P=. ,(b*Q)+{:P
    Q=. Q, <.|(D-*:P)%Q
    whilst. ~:/_2{.P do.
      b=. <.(+/0 _1{P)%{:Q
      P=. P,|(b*{:Q)-{:P
      Q=. Q,|(_2{Q)+b*-/_2{.P
    end.
    f=. y+.x:_2{Q
    if. -. f e. 1,y do. f return. end.
  end.
  1
}}

Task examples:

   task ''
2501: 61 * 41
12851: 71 * 181
13289: 137 * 97
75301: 293 * 257
120787: 43 * 2809
967009: 601 * 1609
997417: 257 * 3881
7091569: 2663 * 2663
13290059: 3119 * 4261
42854447: 9689 * 4423
223553581: 11213 * 19937
2027651281: 46061 * 44021
11111111111: 21649 * 513239
100895598169: 112303 * 898423
1002742628021 was not factored
60012462237239: 6862753 * 8744663
287129523414791: 6059887 * 47381993
9007199254740931: 10624181 * 847801751
11111111111111111: 2071723 * 5363222357
314159265358979323: 317213509 * 990371647
384307168202281507: 415718707 * 924440401
419244183493398773: 48009977 * 8732438749
658812288346769681: 62222119 * 10588072199
922337203685477563: 110075821 * 8379108103
1000000000000000127 was not factored
1152921505680588799: 139001459 * 8294312261
1537228672809128917 was not factored
4611686018427387877 was not factored

where

task=: {{
  for_num. nums do. 
     factor=. x:sqff num
     if. 1=factor do. echo num,&":' was not factored'
     else. echo num,&":': ',factor,&":' * ',":x:num%factor
     end.
  end.
}}

nums=: ".{{)n
        2501
        12851
        13289
        75301
        120787
        967009
        997417
        7091569
        13290059
        42854447
        223553581
        2027651281
        11111111111
        100895598169
        1002742628021
        60012462237239
        287129523414791
        9007199254740931
        11111111111111111
        314159265358979323
        384307168202281507
        419244183493398773
        658812288346769681
        922337203685477563
        1000000000000000127
        1152921505680588799
        1537228672809128917
        4611686018427387877x
}}-.LF

Java

import java.math.BigInteger;
import java.util.List;
import java.util.concurrent.ThreadLocalRandom;

public final class SquareFormFactorization {

	public static void main(String[] args) {		
		ThreadLocalRandom random = ThreadLocalRandom.current();
		final long lowerLimit = 10_000_000_000_000_000L;
		final List<Long> tests = random.longs(20, lowerLimit, 10 * lowerLimit).boxed().toList();

	    for ( long test : tests ) {
	        long factor = squfof(test);	       

	        if ( factor == 0 ) {
	        	System.out.println(test + " - failed to factorise");
	        } else if ( factor == 1 ) {
	        	System.out.println(test + " is a prime number");
	        } else {
	        	System.out.println(test + " = " + factor + " * " + test / factor);
	        }
	        System.out.println();
	   }
	}
	
	private static long squfof(long number) {
		if ( BigInteger.valueOf(number).isProbablePrime(15) ) {
	    	return 1; // Prime number
	    }
		
	    final int sqrt = (int) Math.sqrt(number);
	    if ( sqrt * sqrt == number ) {
	    	return sqrt;
	    }	  
	        
        testValue = number;
        sqrtTestValue = (long) Math.sqrt(testValue);

        // Principal form
        BQF form = new BQF(0, sqrtTestValue, 1);
        form = form.rhoInverse();

		// Search principal cycle
		for ( int i = 0; i < 4 * (long) Math.sqrt(2 * sqrtTestValue); i += 2 ) {
			// Even step
			form = form.rho();				

			long sqrtC = (long) Math.sqrt(form.c);
			if ( sqrtC * sqrtC == form.c ) { // Square form found
				// Inverse square root
				BQF formInverse = new BQF(0, -form.b, sqrtC);
				formInverse = formInverse.rhoInverse();

				// Search ambiguous cycle
				long previousB = 0;
				do {
					previousB = formInverse.b;
					formInverse = formInverse.rho();
				} while ( formInverse.b != previousB );
				
				// Symmetry point
				final long gcd = gcd(number, formInverse.a);
				if ( gcd != 1 ) {
					return gcd;
				}
			}
			
			// Odd step
			form = form.rho();				
		}
		
		if ( number % 2 == 0 ) {
			return 2;
		}	    
		return 0; // Failed to factorise
	}

	private static long gcd(long a, long b) {
		while ( b != 0 ) {
			long temp = a; a = b; b = temp % b;
		}
		return a;
	}	
	
	private static class BQF { // Binary quadratic form
		
		public BQF(long aA, long aB, long aC) {
			a = aA; b = aB; c = aC;
			q = ( sqrtTestValue + b ) / c;
			bb = q * c - b;
		}
		
		public BQF rho() {
			return new BQF(c, bb, a  + q * ( b - bb ));		
		}
		
		public BQF rhoInverse() {
			return new BQF(c, bb, ( testValue - bb * bb ) / c);
		}
		
		private long a, b, c;
		private long q, bb;
		
	}
	
	private static long testValue, sqrtTestValue;

}
Output:
20096060843736547 = 433 * 46411225967059

24628423963378844 = 7 * 3518346280482692

68276045265502398 = 37 * 1845298520689254

61072103663732497 = 8477 * 7204447760261

63462639942509072 = 16 * 3966414996406817

60313009405143787 = 89288189 * 675486983

76093594148871700 = 377900 * 201359074223

31796652636180617 is a prime number

87047981623879461 = 243 * 358222146600327

71567116631895554 = 73 * 980371460710898

50852012325831410 = 2 * 25426006162915705

65816967116185802 = 131280559 * 501345878

89627452852493643 = 31 * 2891208156532053

41735751565855318 = 10004047 * 4171886794

97291513005945602 = 2 * 48645756502972801

88974788272758998 = 59 * 1508047258860322

53903340306287681 = 21727 * 2480938017503

10811459482792395 = 546427 * 19785734385

95115727966103864 = 26105228 * 3643550938

11340988571009785 = 5 * 2268197714201957

jq

Adapted from Wren

Works with gojq, the Go implementation of jq

gojq has support for unbounded-precision integer arithmetic and accordingly the output shown below is from a run thereof; the C implementation of jq produces correct results up to and including [287129523414791,6059887,47381993].

Preliminaries

def gcd(a; b):
  # subfunction expects [a,b] as input
  # i.e. a ~ .[0] and b ~ .[1]
  def rgcd: if .[1] == 0 then .[0]
         else [.[1], .[0] % .[1]] | rgcd
         end;
  [a,b] | rgcd;

# for infinite precision integer-arithmetic
def idivide($p; $q): ($p - ($p % $q)) / $q ;
def idivide($q): (. - (. % $q)) / $q ;

def isqrt:
  def irt:
  . as $x
    | 1 | until(. > $x; . * 4) as $q
    | {$q, $x, r: 0}
    | until( .q <= 1;
        .q |= idivide(4)
        | .t = .x - .r - .q
        | .r |= idivide(2)
        | if .t >= 0
          then .x = .t
          | .r += .q
          else .
          end)
    | .r ;
  if type == "number" and (isinfinite|not) and (isnan|not) and . >= 0
  then irt
  else "isqrt requires a non-negative integer for accuracy" | error
  end ;

The Tasks

def multipliers:
 [
    1, 3, 5, 7, 11, 3*5, 3*7, 3*11, 5*7, 5*11, 7*11, 3*5*7, 3*5*11, 3*7*11, 5*7*11, 3*5*7*11
 ];

# input should be a number
def squfof:
  def toi : floor | tostring | tonumber;
  . as $N
  | (($N|sqrt + 0.5)|toi) as $s
  | if ($s*$s == $N) then $s
    else label $out
    | {}
    | multipliers[] as $multiplier
    | ($N * $multiplier) as $D
        | .P = ($D|isqrt)
        | .Pprev = .P
        | .Pprev as $Po
        | .Qprev = 1
        | .Q = $D - $Po*$Po
        | (($s * 8)|isqrt) as $L
        | (3 * $L) as $B
        | .i = 2
        | .b = 0
        | .q = 0
        | .r = 0
	| .stop = false
        | until( (.i >= $B) or .stop;
            .b = idivide($Po + .P; .Q)
            | .P = .b * .Q - .P
            | .q = .Q
            | .Q = .Qprev + .b * (.Pprev - .P)

            | .r = (((.Q|isqrt) + 0.5)|toi)

            | if ((.i % 2) == 0 and (.r*.r) == .Q) then .stop = true
	      else
                .Qprev = .q
              | .Pprev = .P
              |  .i += 1
	      end )
        | if .i < $B
	  then
            .b = idivide($Po - .P; .r)
	    | .P = .b*.r + .P
            | .Pprev = .P
            | .Qprev = .r
            | .Q = idivide($D - .Pprev*.Pprev; .Qprev)
            | .i = 0
	    | .stop = false
            | until (.stop;
	        .b = idivide($Po + .P; .Q)
                | .Pprev = .P
                | .P = .b * .Q - .P
                | .q = .Q
                | .Q = .Qprev + .b * (.Pprev - .P)
                | .Qprev = .q
                | .i += 1
                | if (.P == .Pprev) then .stop = true else . end )
            | .r = gcd($N; .Qprev)
            | if .r != 1 and .r != $N then .r, break $out else empty end 
          else empty
          end
    end
    // 0 ;

def examples: [
    "2501",
    "12851",
    "13289",
    "75301",
    "120787",
    "967009",
    "997417",
    "7091569",
    "13290059",
    "42854447",
    "223553581",
    "2027651281",
    "11111111111",
    "100895598169",
    "1002742628021",
    "60012462237239",
    "287129523414791",
    "9007199254740931",
    "11111111111111111",
    "314159265358979323",
    "384307168202281507",
    "419244183493398773",
    "658812288346769681",
    "922337203685477563",
    "1000000000000000127",
    "1152921505680588799",
    "1537228672809128917",
    "4611686018427387877"
];
 
"[Integer, Factor, Quotient]"
"---------------------------",
(examples[] as $example
  | ($example|tonumber) as $N
  | ($N | squfof) as $fact
  | if $fact == 0 then "fail"
    else idivide($N; $fact) as $quot
    | [$N, $fact, $quot]
    end
)
Output:
[Integer, Factor, Quotient]
---------------------------
[2501,61,41]
[12851,71,181]
[13289,137,97]
[75301,293,257]
[120787,43,2809]
[967009,1609,601]
[997417,257,3881]
[7091569,2663,2663]
[13290059,3119,4261]
[42854447,9689,4423]
[223553581,11213,19937]
[2027651281,46061,44021]
[11111111111,21649,513239]
[100895598169,112303,898423]
fail
[60012462237239,6862753,8744663]
[287129523414791,6059887,47381993]
[9007199254740931,10624181,847801751]
[11111111111111111,2071723,5363222357]
[314159265358979323,317213509,990371647]
[384307168202281507,415718707,924440401]
[419244183493398773,48009977,8732438749]
[658812288346769681,62222119,10588072199]
[922337203685477563,110075821,8379108103]
[1000000000000000127,111756107,8948056861]
[1152921505680588799,139001459,8294312261]
[1537228672809128917,26675843,57626245319]
[4611686018427387877,343242169,13435662733]

Julia

Modified from Wikipedia's article at [[3]]

function square_form_factor(n::T)::T where T <: Integer
    multiplier = T.([1, 3, 5, 7, 11, 3*5, 3*7, 3*11, 5*7, 5*11, 7*11, 3*5*7, 3*5*11, 3*7*11, 5*7*11, 3*5*7*11])
    s = T(round(sqrt(n)))
    s * s == n && return s
    for k in multiplier
        T != BigInt && n > typemax(T) ÷ k && break
        d = k * n
        p0 = pprev = p = isqrt(d)
        qprev = one(T)
        Q = d - p0 * p0
        l = T(floor(2 * sqrt(2 * s)))
        B, i = 3 * l, 2
        while i < B
            b = (p0 + p) ÷ Q
            p = b * Q - p
            q = Q
            Q = qprev + b * (pprev - p)
            r = T(round(sqrt(Q)))
            iseven(i) && r * r == Q && break
            qprev, pprev = q, p
            i += 1
        end
        i >= B && continue
        b = (p0 - p) ÷ r
        pprev = p = b * r + p
        qprev = r
        Q = (d - pprev * pprev) ÷ qprev
        i = 0
        while true
            b = (p0 + p) ÷ Q
            pprev = p
            p = b * Q - p
            q = Q
            Q = qprev + b * (pprev - p)
            qprev = q
            i += 1
            p == pprev && break
        end
        r = gcd(n, qprev)
        r != 1 && r != n && return r
    end
    return zero(T)
end

println("Integer              Factor   Quotient\n", "-"^45)
@time for n in Int128.([
    2501, 12851, 13289, 75301, 120787, 967009, 997417, 7091569, 13290059, 42854447, 223553581,
    2027651281, 11111111111, 100895598169, 1002742628021, 60012462237239, 287129523414791,
    9007199254740931, 11111111111111111, 314159265358979323, 384307168202281507, 419244183493398773,
    658812288346769681, 922337203685477563, 1000000000000000127, 1152921505680588799,
    1537228672809128917, 4611686018427387877])
    print(rpad(n, 22))
    factr = square_form_factor(n)
    print(rpad(factr, 10))
    println(factr == 0 ? "fail" : n ÷ factr)
end
Output:
Integer              Factor   Quotient
---------------------------------------------
2501                  61        41
12851                 71        181
13289                 137       97
75301                 293       257
120787                43        2809
967009                1609      601
997417                257       3881
7091569               2663      2663
13290059              3119      4261
42854447              9689      4423
223553581             11213     19937
2027651281            46061     44021
11111111111           21649     513239
100895598169          112303    898423
1002742628021         0         fail
60012462237239        6862753   8744663
287129523414791       6059887   47381993
9007199254740931      10624181  847801751
11111111111111111     2071723   5363222357
314159265358979323    317213509 990371647
384307168202281507    415718707 924440401
419244183493398773    48009977  8732438749
658812288346769681    62222119  10588072199
922337203685477563    110075821 8379108103
1000000000000000127   111756107 8948056861
1152921505680588799   139001459 8294312261
1537228672809128917   26675843  57626245319
4611686018427387877   343242169 13435662733
  0.039027 seconds (698 allocations: 38.312 KiB)

Nim

Translation of: Phix
import math, strformat

const M = [uint64 1, 3, 5, 7, 11]

template isqrt(n: uint64): uint64 = uint64(sqrt(float(n)))
template isEven(n: uint64): bool = (n and 1) == 0

proc squfof(n: uint64): uint64 =

  if n.isEven: return 2
  var h = uint64(sqrt(float(n)) + 0.5)
  if h * h == n: return h

  for m in M:
    if m > 1 and (n mod m == 0): return m
    # Check overflow m * n.
    if n > uint64.high div m: break
    let mn = m * n
    var r = isqrt(mn)
    if r * r > mn: dec r
    let rn = r

    # Principal form.
    var b = r
    var a = 1u64
    h = (rn + b) div a * a - b
    var c = (mn - h * h) div a

    for i in 2..<(4 * isqrt(2 * r)):
      # Search principal cycle.
      swap a, c
      var q = (rn + b) div a
      let t = b
      b = q * a - b
      c += q * (t - b)

      if i.isEven:
        r = uint64(sqrt(float(c)) + 0.5)
        if r * r == c: # Square form found?

          # Inverse square root.
          q = (rn - b) div r
          var v = q * r + b
          var w = (mn - v * v) div r

          # Search ambiguous cycle.
          var u = r
          while true:
            swap w, u
            r = v
            q = (rn + v) div u
            v = q * u - v
            if v == r: break
            w += q * (r - v)

          # Symmetry point.
          h = gcd(n, u)
          if h != 1: return h

  result = 1

const Data = [2501u64,
              12851u64,
              13289u64,
              75301u64,
              120787u64,
              967009u64,
              997417u64,
              7091569u64,
              13290059u64,
              42854447u64,
              223553581u64,
              2027651281u64,
              11111111111u64,
              100895598169u64,
              1002742628021u64,
              60012462237239u64,
              287129523414791u64,
              9007199254740931u64,
              11111111111111111u64,
              314159265358979323u64,
              384307168202281507u64,
              419244183493398773u64,
              658812288346769681u64,
              922337203685477563u64,
              1000000000000000127u64,
              1152921505680588799u64,
              1537228672809128917u64,
              4611686018427387877u64]

echo "N                      f          N/f"
echo "======================================"
for n in Data:
  let f = squfof(n)
  let res = if f == 1: "fail" else: &"{f:<10} {n div f}"
  echo &"{n:<22} {res}"
Output:
N                      f          N/f
======================================
2501                   61         41
12851                  71         181
13289                  97         137
75301                  293        257
120787                 43         2809
967009                 1609       601
997417                 257        3881
7091569                2663       2663
13290059               3119       4261
42854447               4423       9689
223553581              11213      19937
2027651281             44021      46061
11111111111            21649      513239
100895598169           112303     898423
1002742628021          fail
60012462237239         6862753    8744663
287129523414791        6059887    47381993
9007199254740931       10624181   847801751
11111111111111111      2071723    5363222357
314159265358979323     317213509  990371647
384307168202281507     415718707  924440401
419244183493398773     48009977   8732438749
658812288346769681     62222119   10588072199
922337203685477563     110075821  8379108103
1000000000000000127    111756107  8948056861
1152921505680588799    139001459  8294312261
1537228672809128917    26675843   57626245319
4611686018427387877    343242169  13435662733

Pascal

Works with: Free Pascal

A console program written in Free Pascal. The code is based on: Jason E. Gower and Samuel S. Wagstaff, jr., "Square form factorization", Mathematics of Computation Volume 77, Number 261, January 2008, Pages 551–588 S 0025-5718(07)02010-8 Article electronically published on May 14, 2007

I'm not sure whether this is the same as reference [1] in the task description; the words "a detailed analysis of SquFoF" appear in the abstract.

The Pascal program includes some small changes to the Gower-Wagstaff algorithm, as noted in the comments. The output shows the successful multiplier (if any) and whether the factor was found in the main or preliminary part of the algorithm.

Further to the advice (on the Discussion page) not to use the Wikipedia version of the algorithm, I tested the Gower-Wagstaff version for all odd composite numbers less than 10^9, and it found a factor in every case. The Wikipedia version failed in 790 cases.

program SquFoF_console;

{$mode objfpc}{$H+}

uses SquFoF_utils;

type TResultKind =
  (rkPrelim, // a factor was found by the preliminary routine
   rkMain,   // a factor was found by the main algorithm
   rkFail);  // no factor was found
type TAlgoResult = record
  Kind : TResultKind;
  Mult : word;
  Factor : UInt64;
end;

// Preliminary to G-W algorithm. Returns D and S of the algorithm.
// Also returns a non-trivial factor if found, else returns factor = 1.
procedure GWPrelim( N : UInt64; // number to be factorized
                    m : word;   // multiplier
                    out D, S, factor : UInt64);
var
  sqflag : boolean;
begin
  D := m*N;
  sqflag := SquFoF_utils.IsSquare( D, S);
  if m = 1 then
    if sqflag then factor := S
              else factor := 1
  else
    factor := GCD( N,m);
end;

// Tries to factorize N by applying Gower-Wagstaff alsorithm to m*N.
function GW_with_multiplier( N : UInt64;
                             m : word) : TAlgoResult;
var
  D, S, P, P_prev, Q, L, B: Uint64;
  r : UInt64;
  i, j, k : integer;
  f, g : UInt64;
  sqrtD : double;
  endCycle : boolean;
// Queue is not much used, so make it a simple array.
type TQueueItem = record
  Left, Right : UInt64;
end;
const QUEUE_CAPACITY = 50; // as suggested by Gower & Wagstaff
var
  queue : array [0..QUEUE_CAPACITY - 1] of TQueueItem;
  queueCount : integer;
begin
  result.Mult := m;

// Filter out special cases (differs from Gower & Wagstaff). Note:
// (1) multiplier m is assumed to be squarefree;
// (2) if we proceed to the main algorithm, mN must not be square
//     (otherwise Q = 0 and division by Q causes an error).
  GWPrelim( N, m, {out} D, S, f);
  if f > 1 then begin
    result.Kind := rkPrelim;
    result.Factor := f;
    exit;
  end;
// Not special, proceed to main algorithm
  result.Kind := rkMain;
  result.Mult := m;
  result.Factor := 1;
  queueCount := 0; // Clear queue
  P := S;
  Q := 1;
  i := -1; // keep i same as in G & W; algo fails if i > B
  sqrtD := SquFoF_utils.FSqrt( D);
//  L := Trunc( 2.0*Sqrt( 2.0*sqrtD));
  L := Trunc( Sqrt( 2.0*sqrtD)); // as in Section 5.2 of G&W paper
  B := 2*L;

  // Start forward cycle
  endCycle := false;
  while not endCycle do begin
    // We update Q here, P at the end of the loop
    Q := (D - P*P) div Q;
    if (not Odd(i)) and SquFoF_utils.IsSquare( Q, r) then begin
      // Q is square for even i.
      // Possibly (?probably?) ends the forward cycle,
      //   but we need to inspect the queue first.
      endCycle := true;
      j := queueCount; // working backwards down the queue
      if r = 1 then begin // the method may fail
        while (j > 0) and (result.Kind <> rkFail) do begin
          dec( j);
          if queue[j].Left = 1 then result.Kind := rkFail;
        end;
        if result.Kind = rkFail then exit;
      end
      else begin // if r > 1
        while (j > 0) and (endCycle) do begin
          dec( j);
          if (queue[j].Left = r)
          and ((P - queue[j].Right) mod r = 0) then begin
              // Deleting up to the *last* item in the list that
              // satisfies the condition, Should it be the *first* ?
              // Delete queue items 0..j inclusive
            inc(j);  k := 0;
            while j < queueCount do begin
              queue[k] := queue[j];
              inc(j);  inc(k);
            end;
            queueCount := k;
            endCycle := false;
          end; // if
        end;
      end;
    end; // if i even and Q square
    if not endCycle then begin
      g := Q div SquFoF_utils.GCD( Q, 2*m);
      if g <= L then begin
        if queueCount < QUEUE_CAPACITY then begin
          with queue[queueCount] do begin
            Left := g;  Right := P mod g;
          end;
          inc( queueCount);
        end
        else begin // queue overflow, fail
          result.Kind := rkFail;
          exit;
        end;
      end;
      inc(i);
      if i > B then begin
        result.Kind := rkFail;
        exit;
      end;
    end;
    P := S - ((S + P) mod Q);
  end; // while not endCycle
  Assert( (D - P*P) mod r = 0); // optional check
  P := S - ((S + P) mod r);
  Q := r;
  // Start backward cycle
  endCycle := false;
  while not endCycle do begin
    P_prev := P;
    Q := (D - P*P) div Q;
    P := S - ((S + P) mod q);
    endCycle := (P = P_prev);
  end; // while not endCycle
  // Finished
  result.Factor := Q div SquFoF_utils.GCD( Q, 2*m);
end;

const NR_RC_VALUES = 28;
RC_VALUES : array [0..NR_RC_VALUES - 1] of UInt64 =
( 2501, 12851, 13289, 75301, 120787, 967009, 997417, 7091569, 13290059,
  42854447, 223553581, 2027651281, 11111111111, 100895598169, 1002742628021,
  60012462237239, 287129523414791, 9007199254740931, 11111111111111111,
  314159265358979323, 384307168202281507, 419244183493398773,
  658812288346769681, 922337203685477563, 1000000000000000127,
  1152921505680588799, 1537228672809128917, 4611686018427387877);

type TMultAndMaxN = record
  Mult : word;   // small multiplier
  MaxN : UInt64; // maximum N for that multiplier (N*multiplier < 2^64)
end;
const NR_MULTIPLIERS = 16;
const MULTIPLIERS : array [0..NR_MULTIPLIERS - 1] of TMultAndMaxN =
((Mult:    1; MaxN: 18446744073709551615),
 (Mult:    3; MaxN:  6148914691236517205),
 (Mult:    5; MaxN:  3689348814741910323),
 (Mult:    7; MaxN:  2635249153387078802),
 (Mult:   11; MaxN:  1676976733973595601),
 (Mult:   15; MaxN:  1229782938247303441),
 (Mult:   21; MaxN:   878416384462359600),
 (Mult:   33; MaxN:   558992244657865200),
 (Mult:   35; MaxN:   527049830677415760),
 (Mult:   55; MaxN:   335395346794719120),
 (Mult:   77; MaxN:   239568104853370800),
 (Mult:  105; MaxN:   175683276892471920),
 (Mult:  165; MaxN:   111798448931573040),
 (Mult:  231; MaxN:    79856034951123600),
 (Mult:  385; MaxN:    47913620970674160),
 (Mult: 1155; MaxN:    15971206990224720));

function GowerWagstaff( N : UInt64) : TAlgoResult;
var
 j : integer;
begin
 j := 0;
 result.Kind := rkFail;
 while (result.Kind = rkFail)
   and (j < NR_MULTIPLIERS)
   and (N <= MULTIPLIERS[j].MaxN) do
 begin
   result := GW_with_multiplier( N, MULTIPLIERS[j].Mult);
   if result.Kind = rkFail then inc(j);
 end;
end;

// Main program
var
  j : integer;
  ar : TAlgoResult;
  kindStr : string;
  N, cofactor : UInt64;
begin
  WriteLn( '              Number Mult  M/P    Factorization');
  for j := 0 to NR_RC_VALUES - 1 do begin
    N := RC_VALUES[j];
    ar := GowerWagstaff( N);
    if ar.Kind = rkFail then
      WriteLn( N:20, '  No factor found')
    else begin
      case ar.Kind of
        rkPrelim: kindStr := 'Prelim';
        rkMain : kindStr := 'Main  ';
      end;
      cofactor := N div ar.Factor;
      Assert( cofactor * ar.Factor = N); // check that all has gone well
      WriteLn( N:20, ar.Mult:5, '  ',
               kindStr:6, ' ', ar.Factor, ' * ', cofactor);
    end;
  end;
end.

unit SquFoF_utils;

{$mode objfpc}{$H+}

interface

// Returns floating-point square root of 64-bit unsigned integer.
function FSqrt( x : UInt64) : double;

// Returns whether a 64-bit unsigned integer is a perfect square.
// In either case, returns floor(sqrt(x)) in the out parameter.
function IsSquare( x : UInt64; out iroot : UInt64) : boolean;

// Returns g.c.d. of 64-bit and 16-bit unsigned integer.
function GCD( u : UInt64; x : word) : word;

implementation

function FSqrt( x : UInt64) : double;
// Both Free Pascal and Delphi 7 seem unreliable when casting
// a 64-bit integer to floating point. We use a workaround.
type TSplitUint64 = packed record case boolean of
  true:  (All : UInt64);
  false: (Lo, Hi : longword); // longword is 32-bit unsigned
end;
var
  temp : TSplitUInt64;
begin
  temp.All := x;
  result := Sqrt( 1.0*temp.Lo + 4294967296.0*temp.Hi);
end;

// Based on Rosetta Code ISqrt, solution for Modula-2.
// Trunc of the f.p. square root won't do, bacause of rounding errors..
function IsSquare( x : UInt64; out iroot : UInt64) : boolean;
var
  Xdiv4, q, r, s, z : UInt64;
begin
  Xdiv4 := X shr 2;
  q := 1;
  while q <= Xdiv4 do q := q shl 2;
  z := x;
  r := 0;
  repeat
    s := q + r;
    r := r shr 1;
    if z >= s then begin
      z := z - s;
      r := r + q;
    end;
    q := q shr 2;
  until q = 0;
  iroot := r;
  result := (z = 0);
end;

function GCD( u : UInt64; x : word) : word;
var
  y, t : word;
begin
  y := u mod x;
  while y <> 0 do begin
    t := x mod y;
    x := y;
    y := t;
  end;
  result := x;
end;
end.
Output:
              Number Mult  M/P    Factorization
                2501    3  Main   61 * 41
               12851    1  Main   71 * 181
               13289    1  Main   97 * 137
               75301    3  Main   293 * 257
              120787    1  Main   43 * 2809
              967009    7  Main   1609 * 601
              997417    1  Main   257 * 3881
             7091569    1  Prelim 2663 * 2663
            13290059    1  Main   3119 * 4261
            42854447    3  Main   4423 * 9689
           223553581    1  Main   11213 * 19937
          2027651281    1  Main   44021 * 46061
         11111111111    3  Main   21649 * 513239
        100895598169   11  Main   898423 * 112303
       1002742628021  No factor found
      60012462237239    1  Main   6862753 * 8744663
     287129523414791    5  Main   6059887 * 47381993
    9007199254740931    1  Main   10624181 * 847801751
   11111111111111111    1  Main   2071723 * 5363222357
  314159265358979323    1  Main   317213509 * 990371647
  384307168202281507    1  Main   415718707 * 924440401
  419244183493398773    1  Main   48009977 * 8732438749
  658812288346769681    1  Main   62222119 * 10588072199
  922337203685477563    1  Main   110075821 * 8379108103
 1000000000000000127    1  Main   111756107 * 8948056861
 1152921505680588799    3  Main   139001459 * 8294312261
 1537228672809128917    3  Main   26675843 * 57626245319
 4611686018427387877    1  Main   343242169 * 13435662733

Perl

Translation of: Raku
Library: ntheory
use strict;
use warnings;
use feature 'say';
use ntheory <is_prime gcd forcomb vecprod>;

my @multiplier;
my @p = <3 5 7 11>;
forcomb { push @multiplier, vecprod @p[@_] } scalar @p;

sub sff {
   my($N) = shift;
   return 1 if is_prime $N;                    # if n is prime
   return sqrt $N if sqrt($N) == int sqrt $N;  # if n is a perfect square

   for my $k (@multiplier) {
      my $P0 = int sqrt($k*$N);  # P[0]=floor(sqrt(N)
      my $Q0 = 1;                # Q[0]=1
      my $Q  = $k*$N - $P0**2;   # Q[1]=N-P[0]^2 & Q[i]
      my $P1 = $P0;              # P[i-1] = P[0]
      my $Q1 = $Q0;              # Q[i-1] = Q[0]
      my $P  = 0;                # P[i]
      my $Qn = 0;                # $P[$i+1];
      my $b  = 0;                # b[i]

      until (sqrt($Q) == int(sqrt($Q))) {           # until Q[i] is a perfect square
         $b = int( int(sqrt($k*$N) + $P1 ) / $Q);   # floor(floor(sqrt(N+P[i-1])/Q[i])
         $P = $b*$Q - $P1;                          # P[i]=b*Q[i]-P[i-1]
         $Qn = $Q1 + $b*($P1 - $P);                 # Q[i+1]=Q[i-1]+b(P[i-1]-P[i])
         ($Q1, $Q, $P1) = ($Q, $Qn, $P);
      }

      $b  = int( int( sqrt($k*$N)+$P ) / $Q );      # b=floor((floor(sqrt(N)+P[i])/Q[0])
      $P1 = $b*$Q0 - $P;                            # P[i-1]=b*Q[0]-P[i]
      $Q  = ( $k*$N - $P1**2 )/$Q0;                 # Q[1]=(N-P[0]^2)/Q[0] & Q[i]
      $Q1 = $Q0;                                    # Q[i-1] = Q[0]

      while () {
         $b  = int( int(sqrt($k*$N)+$P1  ) / $Q );  # b=floor(floor(sqrt(N)+P[i-1])/Q[i])
         $P  = $b*$Q - $P1;                         # P[i]=b*Q[i]-P[i-1]
         $Qn = $Q1 + $b*($P1 - $P);                 # Q[i+1]=Q[i-1]+b(P[i-1]-P[i])
         last if $P == $P1;                         # until P[i+1]=P[i]
         ($Q1, $Q, $P1) = ($Q, $Qn, $P);
      }
      for (gcd $N, $P) { return $_ if $_ != 1 and $_ != $N }
   }
   return 0
}

for my $data (
   11111, 2501, 12851, 13289, 75301, 120787, 967009, 997417,  4558849,  7091569, 13290059,
   42854447, 223553581, 2027651281, 11111111111, 100895598169, 1002742628021, 60012462237239, 
   287129523414791, 11111111111111111, 384307168202281507, 1000000000000000127, 9007199254740931, 
   922337203685477563, 314159265358979323, 1152921505680588799, 658812288346769681,
   419244183493398773, 1537228672809128917) {
   my $v = sff($data);
   if    ($v == 0) { say 'The number ' . $data . ' is not factored.' }
   elsif ($v == 1) { say 'The number ' . $data . ' is a prime.'      }
   else            { say "$data = " . join ' * ', sort {$a <=> $b} $v, int $data/int($v) }
}
Output:
11111 = 41 * 271
2501 = 41 * 61
12851 = 71 * 181
13289 = 97 * 137
75301 = 257 * 293
120787 = 43 * 2809
967009 = 601 * 1609
997417 = 257 * 3881
4558849 = 383 * 11903
7091569 = 2663 * 2663
13290059 = 3119 * 4261
42854447 = 4423 * 9689
223553581 = 11213 * 19937
2027651281 = 44021 * 46061
11111111111 = 21649 * 513239
100895598169 = 112303 * 898423
The number 1002742628021 is a prime.
60012462237239 = 6862753 * 8744663
287129523414791 = 6059887 * 47381993
11111111111111111 = 2071723 * 5363222357
384307168202281507 = 415718707 * 924440401
1000000000000000127 = 111756107 * 8948056861
9007199254740931 = 10624181 * 847801751
922337203685477563 = 110075821 * 8379108103
314159265358979323 = 317213509 * 990371647
1152921505680588799 = 139001459 * 8294312261
658812288346769681 = 62222119 * 10588072199
419244183493398773 = 48009977 * 8732438749
1537228672809128917 = 26675843 * 57626245319

Phix

Translation of: C – (Classical heuristic - fixes the two incorrectly failing cases of the previous version)
--requires(64) -- (decided to limit 32-bit explicitly instead)
constant MxN = power(2,iff(machine_bits()=32?53:63)),
         m = {1, 3, 5, 7, 11}
 
function squfof(atom N)
    -- square form factorization
 
    integer h, a=0, b, c, u=0, v, w, rN, q, r, t
 
    if remainder(N,2)==0 then return 2 end if
    h = floor(sqrt(N) + 0.5)
    if h*h==N then return h end if
 
    for k=1 to length(m) do
        integer mk = m[k]
        if mk>1 and remainder(N,mk)==0 then return mk end if
        //check overflow m * N
        if N>MxN/mk then exit end if
        atom mN = N*mk
        r = floor(sqrt(mN))
        if r*r>mN then r -= 1 end if
        rN = r
 
        //principal form
        {b,a} = {r,1}
        h = floor((rN+b)/a)*a-b
        c = floor((mN-h*h)/a)
 
        for i=2 to floor(sqrt(2*r)) * 4-1 do
            //search principal cycle
 
            {a,c,t} = {c,a,b}
            q = floor((rN+b)/a)
            b = q*a-b
            c += q*(t-b)
 
            if remainder(i,2)==0 then
                r = floor(sqrt(c)+0.5)
 
                if r*r==c then
                    //square form found
 
                    //inverse square root
                    q = floor((rN-b)/r)
                    v = q*r+b
                    w = floor((mN-v*v)/r)
 
                    //search ambiguous cycle
                    u = r
                    while true do
                        {u,w,r} = {w,u,v}
                        q = floor((rN+v)/u)
                        v = q*u-v
                        if v==r then exit end if
                        w += q*(r-v)
                    end while
 
                    //symmetry point
                    h = gcd(N,u)
                    if h!=1 then return h end if
                end if
            end if
        end for
   end for
   return 1
end function
 
constant tests = {2501, 12851, 13289, 75301