basis/math/primes/pollard-rho-brent/pollard-rho-brent.factor

   1 ! Copyright (C) 2021 Doug Coleman.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: arrays kernel make math math.order math.primes
   4 math.ranges random sorting ;
   5 IN: math.primes.pollard-rho-brent
   6
   7 ! https://comeoncodeon.wordpress.com/2010/09/18/pollard-rho-brent-integer-factorization/
   8 :: (brent-factor) ( n -- factor )
   9     n [1..b) random
  10     n [1..b) random
  11     n [1..b) random :> ( y! c m )
  12     1 1 1 :> ( g! r! q! )
  13     0 :> x!
  14     0 :> ys!
  15     [ g 1 = ] [
  16         y x!
  17         r [
  18             y sq n mod c + n mod y!
  19         ] times
  20         0 :> k!
  21         [ k r < g 1 = and ] [
  22             y ys!
  23             m r k - min [
  24                 y sq n mod c + n mod y!
  25                 x y - abs q * n mod q!
  26             ] times
  27             q n gcd nip g!
  28             k m + k!
  29         ] while
  30         r 2 * r!
  31     ] while
  32     g n = [
  33         [ g 1 > not ] [
  34             ys sq n mod c + n mod ys!
  35             x ys - abs n gcd nip g!
  36         ] while
  37     ] when
  38     g ;
  39
  40 : brent-factor ( n -- factor )
  41     dup even? [ drop 2 ] [ (brent-factor) ] if ;
  42
  43 DEFER: pollard-rho-brent-factors
  44
  45 ! 0) can hang if given a large prime, so check prime? first
  46 ! 1) brent-factor can return a composite number
  47 ! 2) also it can hang if you give it a prime number
  48 ! 3) factors can be found in random order
  49 ! therefore we check these conditions
  50 :: (pollard-rho-brent-factors) ( n! -- )
  51     n brent-factor :> factor
  52     ! 1) check for prime
  53     factor prime? [
  54         factor ,
  55     ] [
  56         factor pollard-rho-brent-factors %
  57     ] if
  58     n factor = [
  59         n factor /i
  60         ! 2) check for prime
  61         dup prime? [ , ] [ (pollard-rho-brent-factors) ] if
  62     ] unless ;
  63
  64 : pollard-rho-brent-factors ( n! -- factors )
  65     dup 1 <= [
  66         drop { }
  67     ] [
  68         dup prime? [
  69             1array
  70         ] [
  71             [ (pollard-rho-brent-factors) ] { } make
  72         ] if
  73     ] if natural-sort ;