32-bit Integer Factoring Routine

[aaustin]

Hello,

I have been working on a version for larger integers, but I figured I'd try optimizing the 32-bit size version in case anyone is interested in this stuff. This is one of my favorite methods and is useful when factoring random numbers larger than 128-bits when you want to remove any 64-bit sized primes (or when you want to factor a 128-bit number completely). I'm pretty sure it should be significantly faster than the (direct search/trial divisor) method. At least, for 16-bit numbers or larger I would guess that it should be faster than the simpler method; yet it is simple enough in itself.

FastFactor32 PROC USES EBX ESI EDI, number_to_factor:DWORD
LOCAL limit:DWORD, counter:DWORD
mov esi, number_to_factor
xor eax, eax
cmp esi, 2
jb f_end ;if number is zero or one, it does not have any factors. Return error code (0).
cmp esi, 16
ja @F ;if number is small, quickly factor it using comparisons. Otherwise, jump ahead.
cmp esi, 2
sete dl
cmp esi, 3
sete cl
cmp esi, 5
sete bl
or dl, cl
cmp esi, 7
sete cl
or dl, bl
cmp esi, 11
sete bl
or dl, cl
cmp esi, 13
sete cl
or dl, bl
or dl, cl
movzx eax, dl
jnz f_end ;number is less than 17 and prime? If so, end now with prime return code (1).
add eax, 2
test esi, 1
setnz dl
movzx edx, dl
add eax, edx
jmp f_end
COMMENT !
if number has a factor of two and is not prime, return 2. Otherwise, the number must be 9 or 15.
In that case, return 3.
!
@@: ;for numbers larger than 16.
mov eax, 2
test esi, 1
jz f_end ;number is even, return (2).
bsr eax, esi ;get highest bit set.
shr eax, 2
inc eax ;add 1 for rounding
mov limit, eax ;limit set is very roughly n^(1/4) rounded up to the nearest power of two.
mov eax, esi
shr eax, 1 ; X[0] = floor(n/2)
mul eax
add eax, 1
adc edx, 0
div esi
xor ecx, ecx
mov counter, 1
mov eax, edx ;X[Y+1] = (1 + X[Y]^2) mod n
loop_1:
mov edi, eax
loop_2:
mul eax
add eax, 1
adc edx, 0
div esi
mov eax, edx
mov ebx, edx ; ebx = temp storage register
xor edx, edx
sub eax, edi ;subtract? We really want absolute value.
sbb edx, 0
xor eax, edx
sub eax, edx ; now eax is the absolute value
mov edx, 0
xchg eax, esi
jnz insert_loop
jmp skip_GCD ; I don't really know if this is possible or not. I put it in as a precaution to avoid a Divide Exception (division by 0).
@@: ;GCD internal loop
mov eax, esi
mov esi, edx
xor edx, edx
insert_loop:
div esi
cmp edx, 2
jae @B
mov eax, esi
mov esi, number_to_factor
test edx, edx
jz f_end ;if GCD is 0, return factor in eax
skip_GCD: ;jump here if no GCD was possible.
inc ecx
mov eax, ebx ;restore eax to X[Y].
cmp ecx, count
jb loop_2
mov ebx, count
xor ecx, ecx
add ebx, ebx
cmp ebx, limit
mov count, ebx
jbe loop_1
mov eax, 1 ;factor not found within limit. Number is prime, return 1.
f_end:
FastFactor32 ENDP

COMMENT !
This function uses Brent's Factorization Method which has complexity of roughly n^(1/4) loops.
Polynomial X[Y+1] = (1 + X[Y]^2) mod n
Returns one of the factors. If it is 1, the number is prime. If the number input was 1 or 0,
thr return value is 0 (error).

[MichaelW]

After I corrected the three instances of "count" to "counter" and added a return instruction it seems to work correctly, but for the input values 0-10000 it generates a divide by zero exception for 21, 61, 7177, 7375, and 8375.

Code Select Expand


; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««
    include \masm32\include\masm32rt.inc

    include \masm32\include\debug.inc
    includelib \masm32\lib\debug.lib

    DBGWIN_EXT_INFO = 1
    DBGWIN_DEBUG_ON = 1

    FastFactor32 PROTO :DWORD

; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««
    .data
    .code
; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««
start:
; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««
    TrapException <offset EH>

    mov   ebx, 0
    .WHILE ebx <= 10000
      ;.IF ebx != 21 && ebx != 61 && ebx != 7177 && ebx != 7375 && ebx != 8375
        invoke FastFactor32, ebx
        push eax
        print ustr$(ebx),9
        pop  eax
        print ustr$(eax),13,10
      ;.ENDIF
      add   ebx, 1
    .ENDW

  EH:
   
    inkey "Press any key to exit..."
    exit

; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««

FastFactor32 PROC USES EBX ESI EDI, number_to_factor:DWORD
    LOCAL limit:DWORD, counter:DWORD

    COMMENT !
    This function uses Brent's Factorization Method which has complexity of
    roughly n^(1/4) loops. Polynomial X[Y+1] = (1 + X[Y]^2) mod n
    Returns one of the factors. If it is 1, the number is prime.
    If the number input was 1 or 0, the return value is 0 (error).!

    mov esi, number_to_factor
    xor eax, eax
    cmp esi, 2
    jb f_end ;if number is zero or one, it does not have any factors. Return error code (0).
    cmp esi, 16
    ja @F ;if number is small, quickly factor it using comparisons. Otherwise, jump ahead.
    cmp esi, 2
    sete dl
    cmp esi, 3
    sete cl
    cmp esi, 5
    sete bl
    or dl, cl
    cmp esi, 7
    sete cl
    or dl, bl
    cmp esi, 11
    sete bl
    or dl, cl
    cmp esi, 13
    sete cl
    or dl, bl
    or dl, cl
    movzx eax, dl
    jnz f_end ;number is less than 17 and prime? If so, end now with prime return code (1).
    add eax, 2
    test esi, 1
    setnz dl
    movzx edx, dl
    add eax, edx
    jmp f_end

    COMMENT !
    if number has a factor of two and is not prime, return 2.
    Otherwise, the number must be 9 or 15. In that case, return 3.
    !

  @@:   ;for numbers larger than 16.
    mov eax, 2
    test esi, 1
    jz f_end ;number is even, return (2).
    bsr eax, esi ;get highest bit set.
    shr eax, 2
    inc eax ;add 1 for rounding
    mov limit, eax ;limit set is very roughly n^(1/4) rounded up to the nearest power of two.
    mov eax, esi
    shr eax, 1 ; X[0] = floor(n/2)
    mul eax
    add eax, 1
    adc edx, 0
    div esi
    xor ecx, ecx
    mov counter, 1
    mov eax, edx ;X[Y+1] = (1 + X[Y]^2) mod n
  loop_1:
    mov edi, eax
  loop_2:
    mul eax
    add eax, 1
    adc edx, 0

    ;--------------------------------------------------------------
    ; (address 00401272) divide by zero for N=21,61,7177,7375,8375
    ;--------------------------------------------------------------
    div esi

    mov eax, edx
    mov ebx, edx ; ebx = temp storage register
    xor edx, edx
    sub eax, edi ;subtract? We really want absolute value.
    sbb edx, 0
    xor eax, edx
    sub eax, edx ; now eax is the absolute value
    mov edx, 0
    xchg eax, esi
    jnz insert_loop
    jmp skip_GCD ; I don't really know if this is possible or not. I put it in as a precaution to avoid a Divide Exception (division by 0).
  @@: ;GCD internal loop
    mov eax, esi
    mov esi, edx
    xor edx, edx
    insert_loop:
    div esi
    cmp edx, 2
    jae @B
    mov eax, esi
    mov esi, number_to_factor
    test edx, edx
    jz f_end ;if GCD is 0, return factor in eax
  skip_GCD: ;jump here if no GCD was possible.
    inc ecx
    mov eax, ebx ;restore eax to X[Y].
    cmp ecx, counter
    jb loop_2
    mov ebx, counter
    xor ecx, ecx
    add ebx, ebx
    cmp ebx, limit
    mov counter, ebx
    jbe loop_1
    mov eax, 1 ;factor not found within limit. Number is prime, return 1.
  f_end:
    ret
FastFactor32 ENDP

; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««
end start


[aaustin]

Hmmm....

Thanks for testing it Micheal. I didn't have a chance to last night.

I didn't notice the count/counter discrepancy, oops. Thanks for correcting it. I usually don't put a return statement in my functions. It's easier to specify a language type (such as STDCALL), that way, after cleaning up the stack it automatically inserts a, "ret 4" instruction and takes care of near or far distance automatically (if you put the ret instruction in front of the end of the function, I don't think the local variables get cleaned up).

I know what's wrong. The restoration of esi ("mov esi , number_to_factor") should occur immediately after the "skip_GCD" jump point. If you make that change, the code should work without exceptions. That was my goof. Thanks for finding it!

[aaustin]

Ok,

After doing tests myself, I see that the numbers you mentioned will create an infinite loop that will cause the numbers to falsely be noted as prime (in the case of 21, 7375 and 8375. 61 and 7177 are correctly noted prime.) I need a way to fix this (most likely a retry using different seed input values). I'll look into a way to fix it.

[MichaelW]

FWIW, after a little more testing I determined that a skip_GCD jump takes place before each of the exceptions (N = 21, 61, 7177, 7375, and 8375) and without triggering an exception for N = 407, 5707, 7943, and 8131. ESI is zero each time the jump is taken, the skip_GCD code does not alter ESI before conditionally jumping to Loop1 or Loop2, neither of which alter ESI ahead of the faulting divide.

[aaustin]

making the following changes will fix the algorithm:
...
    sub eax, edx ; now eax is the absolute value
    xor edx, edx
    cmp eax, 2 ;make sure eax is greater than 2, otherwise the GCD will not work.
    xchg eax, esi ;esi changed
    jae insert_loop
    jmp skip_GCD ; I don't really know if this is possible or not. I put it in as a precaution to avoid a Divide Exception (division by 0).
  @@: ;GCD internal loop
    mov eax, esi
    mov esi, edx
    xor edx, edx
    insert_loop:
    div esi
    cmp edx, 2
    jae @B
    mov eax, esi
    test edx, edx
    jz f_end ;if GCD is 0, return factor in eax
  skip_GCD: ;jump here if no GCD was possible.
    mov esi, number_to_factor ;esi restored properly
    inc ecx
    mov eax, ebx ;restore eax to X[Y].
    cmp ecx, counter
    jb loop_2
...

This should fix the algorithm. The problem with 21 is that X[0] is not compared with X[1] and GCDed. If you correct this along with the last tidbit of code I posted, it may fix the entire procedure.

[MichaelW]

I'm not sure I made the changes correctly. The modified procedure has no problem with divide exceptions for input values from 0 to 1,000,000, so I think that problem is probably fixed. But for input values from 2 to 1000, it correctly identifies 168 primes, and incorrectly identifies 89 composite numbers as prime, starting with 21,25,65,95,115,125,145,155,169...

Code Select Expand


; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««
    include \masm32\include\masm32rt.inc

    include \masm32\include\debug.inc
    includelib \masm32\lib\debug.lib

    DBGWIN_EXT_INFO = 1
    DBGWIN_DEBUG_ON = 1

    FastFactor32 PROTO :DWORD

; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««
    .data     
      primes  dd 0002,0003,0005,0007,0011,0013,0017,0019,0023,0029
              dd 0031,0037,0041,0043,0047,0053,0059,0061,0067,0071
              dd 0073,0079,0083,0089,0097,0101,0103,0107,0109,0113
              dd 0127,0131,0137,0139,0149,0151,0157,0163,0167,0173
              dd 0179,0181,0191,0193,0197,0199,0211,0223,0227,0229
              dd 0233,0239,0241,0251,0257,0263,0269,0271,0277,0281
              dd 0283,0293,0307,0311,0313,0317,0331,0337,0347,0349
              dd 0353,0359,0367,0373,0379,0383,0389,0397,0401,0409
              dd 0419,0421,0431,0433,0439,0443,0449,0457,0461,0463
              dd 0467,0479,0487,0491,0499,0503,0509,0521,0523,0541
              dd 0547,0557,0563,0569,0571,0577,0587,0593,0599,0601
              dd 0607,0613,0617,0619,0631,0641,0643,0647,0653,0659
              dd 0661,0673,0677,0683,0691,0701,0709,0719,0727,0733
              dd 0739,0743,0751,0757,0761,0769,0773,0787,0797,0809
              dd 0811,0821,0823,0827,0829,0839,0853,0857,0859,0863
              dd 0877,0881,0883,0887,0907,0911,0919,0929,0937,0941
              dd 0947,0953,0967,0971,0977,0983,0991,0997,1009,1013
              dd 1019,1021,1031,1033,1039,1049,1051,1061,1063,1069
              dd 1087,1091,1093,1097,1103,1109,1117,1123,1129,1151
              dd 1153,1163,1171,1181,1187,1193,1201,1213,1217,1223
              dd 1229,1231,1237,1249,1259,1277,1279,1283,1289,1291
              dd 1297,1301,1303,1307,1319,1321,1327,1361,1367,1373
              dd 1381,1399,1409,1423,1427,1429,1433,1439,1447,1451
              dd 1453,1459,1471,1481,1483,1487,1489,1493,1499,1511
              dd 1523,1531,1543,1549,1553,1559,1567,1571,1579,1583
              dd 1597,1601,1607,1609,1613,1619,1621,1627,1637,1657
              dd 1663,1667,1669,1693,1697,1699,1709,1721,1723,1733
              dd 1741,1747,1753,1759,1777,1783,1787,1789,1801,1811
              dd 1823,1831,1847,1861,1867,1871,1873,1877,1879,1889
              dd 1901,1907,1913,1931,1933,1949,1951,1973,1979,1987
              dd 1993,1997,1999,2003,2011,2017,2027,2029,2039,2053
              dd 2063,2069,2081,2083,2087,2089,2099,2111,2113,2129
              dd 2131,2137,2141,2143,2153,2161,2179,2203,2207,2213
              dd 2221,2237,2239,2243,2251,2267,2269,2273,2281,2287
              dd 2293,2297,2309,2311,2333,2339,2341,2347,2351,2357
              dd 2371,2377,2381,2383,2389,2393,2399,2411,2417,2423
              dd 2437,2441,2447,2459,2467,2473,2477,2503,2521,2531
              dd 2539,2543,2549,2551,2557,2579,2591,2593,2609,2617
              dd 2621,2633,2647,2657,2659,2663,2671,2677,2683,2687
              dd 2689,2693,2699,2707,2711,2713,2719,2729,2731,2741
              dd 2749,2753,2767,2777,2789,2791,2797,2801,2803,2819
              dd 2833,2837,2843,2851,2857,2861,2879,2887,2897,2903
              dd 2909,2917,2927,2939,2953,2957,2963,2969,2971,2999
              dd 3001,3011,3019,3023,3037,3041,3049,3061,3067,3079
              dd 3083,3089,3109,3119,3121,3137,3163,3167,3169,3181
              dd 3187,3191,3203,3209,3217,3221,3229,3251,3253,3257
              dd 3259,3271,3299,3301,3307,3313,3319,3323,3329,3331
              dd 3343,3347,3359,3361,3371,3373,3389,3391,3407,3413
              dd 3433,3449,3457,3461,3463,3467,3469,3491,3499,3511
              dd 3517,3527,3529,3533,3539,3541,3547,3557,3559,3571
              dd 3581,3583,3593,3607,3613,3617,3623,3631,3637,3643
              dd 3659,3671,3673,3677,3691,3697,3701,3709,3719,3727
              dd 3733,3739,3761,3767,3769,3779,3793,3797,3803,3821
              dd 3823,3833,3847,3851,3853,3863,3877,3881,3889,3907
              dd 3911,3917,3919,3923,3929,3931,3943,3947,3967,3989
              dd 4001,4003,4007,4013,4019,4021,4027,4049,4051,4057
              dd 4073,4079,4091,4093,4099,4111,4127,4129,4133,4139
              dd 4153,4157,4159,4177,4201,4211,4217,4219,4229,4231
              dd 4241,4243,4253,4259,4261,4271,4273,4283,4289,4297
              dd 4327,4337,4339,4349,4357,4363,4373,4391,4397,4409
              dd 4421,4423,4441,4447,4451,4457,4463,4481,4483,4493
              dd 4507,4513,4517,4519,4523,4547,4549,4561,4567,4583
              dd 4591,4597,4603,4621,4637,4639,4643,4649,4651,4657
              dd 4663,4673,4679,4691,4703,4721,4723,4729,4733,4751
              dd 4759,4783,4787,4789,4793,4799,4801,4813,4817,4831
              dd 4861,4871,4877,4889,4903,4909,4919,4931,4933,4937
              dd 4943,4951,4957,4967,4969,4973,4987,4993,4999,5003
              dd 5009,5011,5021,5023,5039,5051,5059,5077,5081,5087
              dd 5099,5101,5107,5113,5119,5147,5153,5167,5171,5179
              dd 5189,5197,5209,5227,5231,5233,5237,5261,5273,5279
              dd 5281,5297,5303,5309,5323,5333,5347,5351,5381,5387
              dd 5393,5399,5407,5413,5417,5419,5431,5437,5441,5443
              dd 5449,5471,5477,5479,5483,5501,5503,5507,5519,5521
              dd 5527,5531,5557,5563,5569,5573,5581,5591,5623,5639
              dd 5641,5647,5651,5653,5657,5659,5669,5683,5689,5693
              dd 5701,5711,5717,5737,5741,5743,5749,5779,5783,5791
              dd 5801,5807,5813,5821,5827,5839,5843,5849,5851,5857
              dd 5861,5867,5869,5879,5881,5897,5903,5923,5927,5939
              dd 5953,5981,5987,6007,6011,6029,6037,6043,6047,6053
              dd 6067,6073,6079,6089,6091,6101,6113,6121,6131,6133
              dd 6143,6151,6163,6173,6197,6199,6203,6211,6217,6221
              dd 6229,6247,6257,6263,6269,6271,6277,6287,6299,6301
              dd 6311,6317,6323,6329,6337,6343,6353,6359,6361,6367
              dd 6373,6379,6389,6397,6421,6427,6449,6451,6469,6473
              dd 6481,6491,6521,6529,6547,6551,6553,6563,6569,6571
              dd 6577,6581,6599,6607,6619,6637,6653,6659,6661,6673
              dd 6679,6689,6691,6701,6703,6709,6719,6733,6737,6761
              dd 6763,6779,6781,6791,6793,6803,6823,6827,6829,6833
              dd 6841,6857,6863,6869,6871,6883,6899,6907,6911,6917
              dd 6947,6949,6959,6961,6967,6971,6977,6983,6991,6997
              dd 7001,7013,7019,7027,7039,7043,7057,7069,7079,7103
              dd 7109,7121,7127,7129,7151,7159,7177,7187,7193,7207
              dd 7211,7213,7219,7229,7237,7243,7247,7253,7283,7297
              dd 7307,7309,7321,7331,7333,7349,7351,7369,7393,7411
              dd 7417,7433,7451,7457,7459,7477,7481,7487,7489,7499
              dd 7507,7517,7523,7529,7537,7541,7547,7549,7559,7561
              dd 7573,7577,7583,7589,7591,7603,7607,7621,7639,7643
              dd 7649,7669,7673,7681,7687,7691,7699,7703,7717,7723
              dd 7727,7741,7753,7757,7759,7789,7793,7817,7823,7829
              dd 7841,7853,7867,7873,7877,7879,7883,7901,7907,7919
      nP      dd 0
      nC      dd 0
    .code
; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««
start:
; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««
    TrapException <offset EH>

    print "Return value for input 0 = "
    invoke FastFactor32, 0
    print ustr$(eax),13,10
    print "Return value for input 1 = "
    invoke FastFactor32, 1   
    print ustr$(eax),13,10

    inkey "Press any key to continue..."

    mov esi, OFFSET primes
    mov ebx, 2
    .WHILE ebx <= 1000
      invoke FastFactor32, ebx
      .IF eax == 1
          print ustr$(ebx)
          xor edi, edi
          .WHILE edi < 1000
            .IF [esi+edi*4] == ebx
              print " P"
              inc nP
              .BREAK
            .ENDIF
            inc edi
          .ENDW
          .IF edi == 1000
            print " C"
            inc nC
          .ENDIF
        print chr$(13,10)
      .ENDIF
      add ebx, 1
    .ENDW
    print "nP = "
    print ustr$(nP),13,10
    print "nC = "
    print ustr$(nC),13,10

  EH:
   
    inkey "Press any key to exit..."
    exit

; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««

FastFactor32 PROC USES EBX ESI EDI, number_to_factor:DWORD
    LOCAL limit:DWORD, counter:DWORD

    COMMENT !
    This function uses Brent's Factorization Method which has complexity of
    roughly n^(1/4) loops. Polynomial X[Y+1] = (1 + X[Y]^2) mod n
    Returns one of the factors. If it is 1, the number is prime.
    If the number input was 1 or 0, the return value is 0 (error).!

    mov esi, number_to_factor
    xor eax, eax
    cmp esi, 2
    jb f_end ;if number is zero or one, it does not have any factors. Return error code (0).
    cmp esi, 16
    ja @F ;if number is small, quickly factor it using comparisons. Otherwise, jump ahead.
    cmp esi, 2
    sete dl
    cmp esi, 3
    sete cl
    cmp esi, 5
    sete bl
    or dl, cl
    cmp esi, 7
    sete cl
    or dl, bl
    cmp esi, 11
    sete bl
    or dl, cl
    cmp esi, 13
    sete cl
    or dl, bl
    or dl, cl
    movzx eax, dl
    jnz f_end ;number is less than 17 and prime? If so, end now with prime return code (1).

    add eax, 2
    test esi, 1
    setnz dl
    movzx edx, dl
    add eax, edx
    jmp f_end

    COMMENT !
    if number has a factor of two and is not prime, return 2.
    Otherwise, the number must be 9 or 15. In that case, return 3.
    !

  @@:   ;for numbers larger than 16.
    mov eax, 2
    test esi, 1
    jz f_end ;number is even, return (2).
    bsr eax, esi ;get highest bit set.
    shr eax, 2
    inc eax ;add 1 for rounding
    mov limit, eax ;limit set is very roughly n^(1/4) rounded up to the nearest power of two.
    mov eax, esi
    shr eax, 1 ; X[0] = floor(n/2)
    mul eax
    add eax, 1
    adc edx, 0
    div esi
    xor ecx, ecx
    mov counter, 1
    mov eax, edx ;X[Y+1] = (1 + X[Y]^2) mod n
  loop_1:
    mov edi, eax
  loop_2:
    mul eax
    add eax, 1
    adc edx, 0
    div esi
    mov eax, edx
    mov ebx, edx ; ebx = temp storage register
    xor edx, edx
    sub eax, edi ;subtract? We really want absolute value.
    sbb edx, 0
    xor eax, edx
  ;------------------------------------------------
    sub eax, edx ; now eax is the absolute value
    xor edx, edx
    cmp eax, 2 ;make sure eax is greater than 2, otherwise the GCD will not work.
    xchg eax, esi ;esi changed
    jae insert_loop
    jmp skip_GCD ; I don't really know if this is possible or not. I put it in as a precaution to avoid a Divide Exception (division by 0).
  @@: ;GCD internal loop
    mov eax, esi
    mov esi, edx
    xor edx, edx
    insert_loop:
    div esi
    cmp edx, 2
    jae @B
    mov eax, esi
    test edx, edx
    jz f_end ;if GCD is 0, return factor in eax
  skip_GCD: ;jump here if no GCD was possible.
    mov esi, number_to_factor ;esi restored properly
    inc ecx
    mov eax, ebx ;restore eax to X[Y].
    cmp ecx, counter
    jb loop_2
  ;------------------------------------------------
    mov ebx, counter
    xor ecx, ecx
    add ebx, ebx
    cmp ebx, limit
    mov counter, ebx
    jbe loop_1
    mov eax, 1 ;factor not found within limit. Number is prime, return 1.
  f_end:
    ret
FastFactor32 ENDP

; «««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««
end start




[attachment deleted by admin]

[aaustin]

well, for 21, if you compare X[0] to X[1] as I suggested:

X[0] = floor(21/2) = 10
X[1] = 1+10^2 mod 21 = 17
17-10 = 7
GCD(21, 7) = 7

Doing this comparison will fix some of the incorrect primes. As for some of the numbers, I've noticed that doing an additional GCD of (X[?], n) will fix a lot of the problems (unfortunately that doubles the work). I don't think just doing that will totally fix the problem, but it will probably reduce false prime detection significantly.

25
X[0] = 12
X[1] = 20
GCD(25, 20) = 5

65
X[0] = 32
X[1] = 50
GCD(65, 50) = 5

95
X[0] = 47
X[1] = 25
GCD(95,25) = 5

115
X[0] = 57
X[1] = 5
GCD(115, 5) = 5

125
X[0] = 62
X[1] = 95
GCD(125, 95) = 5

145
X[0] = 72
X[1] = 110
GCD(145, 110) = 5

155
X[0] = 77
X[1] = 40
GCD(155, 40) = 5

169
X[0] = 84
X[1] = 128
X[2] = 161
X[3] = 65
GCD(169, 65) = 13

[roticv]

I suspect there is something wrong with your Brent's Factorization Method. It looks different from what I have found online...

[aaustin]

Brent's factorization method is a form of pollard rho factorization. If what you've found online does look significantly different, then I suspect that it is not really Brent's factorization method but some other form. Mine is a base 2 algorithm (other bases are allowed).