Square form factorization

A just for fun snail ..

References: Algorithm, C program example from the WP page and also the pseudocode from here.

# 20210325 Raku programming solution

my @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 );

sub circumfix:<โŒŠ โŒ‹>{ $^n.floor }; sub prefix:<โˆš>{ $^n.sqrt }; # just for fun

sub SQUFOF ( \๐‘ ) {  

   return  1 if ๐‘.is-prime;     # if n is prime             return  1
   return โˆš๐‘ if โˆš๐‘ == Int(โˆš๐‘);  # if n is a perfect square  return โˆš๐‘

   for @multiplier -> \๐‘˜ { 
      my \Pโ‚’ = $ = โŒŠ โˆš(๐‘˜*๐‘) โŒ‹;          # P[0]=floor(โˆšN)
      my \Qโ‚’ = $ = 1 ;                  # Q[0]=1
      my \Q = $ =  ๐‘˜*๐‘ - Pโ‚’ยฒ;           # Q[1]=N-P[0]^2 & Q[i]
      my \Pโ‚šแตฃโ‚‘แตฅ = $ = Pโ‚’;               # P[i-1] = P[0]
      my \Qโ‚šแตฃโ‚‘แตฅ = $ = Qโ‚’;               # Q[i-1] = Q[0] 
      my \P = $ = 0;                    # P[i]
      my \Qโ‚™โ‚‘โ‚“โ‚œ = $ = 0;                # P[i+1]
      my \b = $ = 0;                    # b[i]
                                        # i = 1
      repeat until โˆšQ == Int(โˆšQ) {      # until Q[i] is a perfect square
         b = โŒŠโŒŠ โˆš(๐‘˜*๐‘) + Pโ‚šแตฃโ‚‘แตฅ โŒ‹ / Q โŒ‹;    # floor(floor(โˆšN+P[i-1])/Q[i])
         P = b*Q - Pโ‚šแตฃโ‚‘แตฅ;                  # P[i]=b*Q[i]-P[i-1] 
         Qโ‚™โ‚‘โ‚“โ‚œ = Qโ‚šแตฃโ‚‘แตฅ + b*(Pโ‚šแตฃโ‚‘แตฅ - P);    # Q[i+1]=Q[i-1]+b(P[i-1]-P[i]) 
         ( Qโ‚šแตฃโ‚‘แตฅ,  Q, Pโ‚šแตฃโ‚‘แตฅ ) = Q, Qโ‚™โ‚‘โ‚“โ‚œ, P;  # i++
      } 

      b     = โŒŠ โŒŠ โˆš(๐‘˜*๐‘)+P โŒ‹ / Q โŒ‹;     # b=floor((floor(โˆšN)+P[i])/Q[0])
      Pโ‚šแตฃโ‚‘แตฅ = b*Qโ‚’ - P;                 # P[i-1]=b*Q[0]-P[i]
      Q     = ( ๐‘˜*๐‘ - Pโ‚šแตฃโ‚‘แตฅยฒ )/Qโ‚’;      # Q[1]=(N-P[0]^2)/Q[0] & Q[i]
      Qโ‚šแตฃโ‚‘แตฅ = Qโ‚’;                       # Q[i-1] = Q[0]
                                        # i = 1
      loop {                            # repeat
         b  = โŒŠ โŒŠ โˆš(๐‘˜*๐‘)+Pโ‚šแตฃโ‚‘แตฅ โŒ‹ / Q โŒ‹;    # b=floor(floor(โˆšN)+P[i-1])/Q[i])
         P  = b*Q - Pโ‚šแตฃโ‚‘แตฅ;                 # P[i]=b*Q[i]-P[i-1]
         Qโ‚™โ‚‘โ‚“โ‚œ = Qโ‚šแตฃโ‚‘แตฅ + b*(Pโ‚šแตฃโ‚‘แตฅ - P);    # Q[i+1]=Q[i-1]+b(P[i-1]-P[i])
         last if (P == Pโ‚šแตฃโ‚‘แตฅ);          # until P[i+1]=P[i]
         ( Qโ‚šแตฃโ‚‘แตฅ,  Q, Pโ‚šแตฃโ‚‘แตฅ ) = Q, Qโ‚™โ‚‘โ‚“โ‚œ, P; # i++ 
      }  
      given ๐‘ gcd P { return $_ if $_ != 1|๐‘ }    
   }  # gcd(N,P[i]) (if != 1 or N) is a factor of N, otherwise try next k
   return 0 # give up
}

race for ( 
   11111, # wikipedia.org/wiki/Shanks%27s_square_forms_factorization#Example
   4558849, # example from talk page
   # all of the rest are taken from the FreeBASIC entry
   2501,12851,13289,75301,120787,967009,997417,7091569,13290059,
   42854447,223553581,2027651281,11111111111,100895598169,1002742628021,
   # time hoarders 
   60012462237239, # = 6862753 * 8744663                     15s 
   287129523414791, # = 6059887 * 47381993                   80s
   11111111111111111, # = 2071723 * 5363222357                2m
   384307168202281507, # = 415718707 * 924440401              5m
   1000000000000000127, # = 111756107 * 8948056861           12m
   9007199254740931, # = 10624181 * 847801751                17m
   922337203685477563, # = 110075821 * 8379108103            41m
   314159265358979323, # = 317213509 * 990371647             61m
   1152921505680588799, # = 139001459 * 8294312261           93m
   658812288346769681, # = 62222119 * 10588072199           112m
   419244183493398773, # = 48009977 * 8732438749            135m
   1537228672809128917, # = 26675843 * 57626245319          254m
   # don't know how to handle this one
   # for 1e-323, 1e-324 { my $*TOLERANCE = $_ ; 
   #     say 4611686018427387877.sqrt โ‰… 4611686018427387877.sqrt.Int }
   # skip the perfect square check and start k with 3 to get the following 
   # 4611686018427387877, # = 343242169 * 13435662733       217m
) -> \data {
   given data.&SQUFOF {
      when 0  { say "The number ", data, " is not factored." }
      when 1  { say "The number ", data, " is a prime." }
      default { say data, " = ", $_, " * ", data div $_.Int }
   }
}

Output:

11111 = 41 * 271 
4558849 = 383 * 11903 
2501 = 61 * 41 
12851 = 71 * 181 
13289 = 97 * 137 
75301 = 293 * 257 
120787 = 43 * 2809 
967009 = 601 * 1609 
997417 = 257 * 3881 
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 

Use NativeCall

Use idea similar to the second approach from here, by compiling the C (Classical heuristic) entry to a shared library and then make use of the squfof routine.

squfof.raku

# 20210326 Raku programming solution

use NativeCall;

constant LIBSQUFOF = '/home/user/LibSQUFOF.so';

sub squfof(uint64 $n) returns uint64 is native(LIBSQUFOF) { * };

race for (
   11111, # wikipedia.org/wiki/Shanks%27s_square_forms_factorization#Example
   4558849, # example from talk page
   # all of the rest are taken from the FreeBASIC entry
   2501,12851,13289,75301,120787,967009,997417,7091569,13290059,
   42854447,223553581,2027651281,11111111111,100895598169,1002742628021,
   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
   4611686018427387877, # = 343242169 * 13435662733
) -> \data {
   given squfof(data) { say data, " = ", $_, " * ", data div $_ }
}

Output:

gcc -Wall -fPIC -shared -o LibSQUFOF.so squfof.c
<p>file LibSQUFOF.so
LibSQUFOF.so: ELF 64-bit LSB shared object, x86-64, version 1 (SYSV), dynamically linked, BuildID[sha1]=f7f9864e5508ba1de80cbc6e2f4d789f4ec448e9, not stripped
raku -c squfof.raku && time ./squfof.raku
Syntax OK
11111 = 41 * 271
4558849 = 383 * 11903
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 = 1 * 1002742628021
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
4611686018427387877 = 343242169 * 13435662733
</p><p>real    0m0.784s
user    0m0.716s
sys     0m0.056s
</p>

Last updated

Was this helpful?