Yafu

[bsquared]

... in the usual place.

This version has some cool new stuff. bchaffin contributed code to run gmp-ecm over multiple threads - now one can do N times as much ECM pretesting in the same amount of time when running multi-threaded. WraithX contributed code to grab the input expression and use it in his special output files rather than the expanded decimal form (e.g. 2^32-1 would be put in log files rather than 4294967295). Brian Gladman contributed assembler files to support 64 bit mod operations in x64 builds. Lastly, I was able to get mingw64 builds to work, which together with a largish pile of new assembly code results in much faster SIQS code for 64 bit platforms - up to 20% in some cases. Linux people should see similar improvements. My 2.93 GHz i7 can now factor my test C85 in 7.5 minutes, single threaded

Additional details:
* the parallel gmp-ecm code relies on shared memory and fork(). Unfortunately, that means that this is only available for *nix platforms.
* A big piece of the SIQS speed improvement comes from resieving of "medium" sized factor base primes. A big thanks to fivemack and R.D. Silverman for helping me to finally understand what it was all about (towards the end of this thread)

Cheers,
- ben.

[alexhiggins732]

Not sure if this has been fixed in the latest version or not. Using 1.22 there is a bug with the factorization of certain numbers.

To reproduce, take N = the square of a Mersenne number larger than [sqrt(MP)] mod MP. Eg, N=(2^536-1)^2 mod (2^1061-1), or N=(2^537-1)^2 mod (2^1061-1), or N=(2^538-1)^2 mod (2^1061-1)... etc.

The modulus will be half of a perfect square.

When trying to factor such N, or even 2N (which would be a perfect square) YAFU Displays the message:

Composite result found, starting re-factorization
factoring

The factors found before the re-factorization are displayed correct but factors after the factorization is restarted are displayed only once. (Since this is half of a square there should be an even number of each factor).

Here is an example:

>> factor( (((2^536)-1)^2) % ((2^1061)-1))
factoring 2470730631192756571685734212877408533319783322316187968223893530608280
51230463069936475077760543364862288913408589858290270762618879142427816178466724
53431386454091076181222520467347973217390137018275366108375849942122621867807049
86963140774948747879449495438661420404801673551395242809808659633830691291835241
2782102528
div: primes less than 10000
fmt: 1000000 iterations
Composite result found, starting re-factorization
factoring 4714479403884569277363099893868183707610867700925548223796814958268782
068604900555166186853844841268392966736734069296073513841042481
div: primes less than 10000
fmt: 1000000 iterations
rho: x^2 + 1, starting 1000 iterations on C133
rho: x^2 + 1, starting 1000 iterations on C128
rho: x^2 + 1, starting 1000 iterations on C123
rho: x^2 + 1, starting 1000 iterations on C118
rho: x^2 + 1, starting 1000 iterations on C113
rho: x^2 + 1, starting 1000 iterations on C108
rho: x^2 + 1, starting 1000 iterations on C104
rho: x^2 + 1, starting 1000 iterations on C97
rho: x^2 + 1, starting 1000 iterations on C89
rho: x^2 + 3, starting 1000 iterations on C89
rho: x^2 + 3, starting 1000 iterations on C81
rho: x^2 + 2, starting 1000 iterations on C81
pp1: starting B1 = 20K, B2 = gmp-ecm default on C81
pp1: starting B1 = 20K, B2 = gmp-ecm default on C81
pp1: starting B1 = 20K, B2 = gmp-ecm default on C81
pm1: starting B1 = 100K, B2 = gmp-ecm default on C81
ecm: 25 curves on C81 input, at B1 = 2K, B2 = gmp-ecm default
ecm: 5 curves on C81 input, at B1 = 11K, B2 = gmp-ecm default
Total factoring time = 5.6719 seconds
Total factoring time = 5.8281 seconds

***factors found***
P1 = 2
P1 = 2
P1 = 2
P1 = 2
P1 = 2
P1 = 2
P1 = 2
P1 = 2
P1 = 2
P1 = 2
P1 = 2
P1 = 7
P1 = 7
P1 = 7
P1 = 7
P2 = 31
P2 = 31
P2 = 71
P2 = 71
P3 = 127
P3 = 127
P3 = 151
P3 = 151
P3 = 337
P3 = 337
P3 = 601
P3 = 601
P4 = 1801
P4 = 1801
P4 = 4201
P4 = 4201
P4 = 7351
P4 = 7351
PRP6 = 106681 <-- Factors below this are only shown once.
PRP6 = 100801
PRP6 = 122921
PRP6 = 152041
PRP5 = 39551
PRP5 = 29191
PRP8 = 10567201
PRP9 = 181165951
PRP8 = 60816001
PRP24 = 535347624791488552837151
PRP57 = 325985508875527587669607097222667557116221139090131514801
ans = 1

Here is the correct factorization, using Dario's ecm applet and inserting the factors --> www.alpertron.com.ar/ECM.HTM

24 707306 311927 565716 857342 128774 085333 197833 223161 879682 238935
306082 805123 046306 993647 507776 054336 486228 891340 858985 829027 076261
887914 242781 617846 672453 431386 454091 076181 222520 467347 973217 390137
018275 366108 375849 942122 621867 807049 869631 407749 487478 794494 954386
614204 048016 735513 952428 098086 596338 306912 918352 412782 102528
= 2 ^ 11 x
7 ^ 4 x
31 ^ 2 x
71 ^ 2 x
127 ^ 2 x
151 ^ 2 x
337 ^ 2 x
601 ^ 2 x
1801 ^ 2 x
4201 ^ 2 x
7351 ^ 2 x
29191 ^ 2 x
39551 ^ 2 x
100801 ^ 2 x
106681 ^ 2 x
122921 ^ 2 x
152041 ^ 2 x
10 567201 ^ 2 x
60 816001 ^ 2 x
181 165951 ^ 2 x
535347 624791 488552 837151 ^ 2 x
325 985508 875527 587669 607097 222667 557116 221139 090131 514801 ^ 2
Number of divisors: 209207 064060
Sum of divisors: 61 635649 470645 320448 927154 220340 253567 976068 388115
429802 242295 103264 789782 628511 075563 142521 434123 144593 673113 150774
244626 607968 427317 961843 751365 536296 822823 875775 212197 930336 840665
903721 912257 427209 329934 656026 114668 900455 551671 596235 435297 212026
147658 993655 578412 220234 033377 903844 302842 158009 995466 594968 300310
968645
Euler's Totient: 9 900811 595773 482739 754491 145945 316486 735981 308425
186148 164446 700676 148732 384132 330414 020041 175142 590822 163841 695079
301611 325028 456791 100843 561490 412864 403721 057504 726139 431010 496015
191329 380914 605509 068662 571201 016317 347822 168715 195190 726637 971205
049789 158287 344030 882177 382324 562129 059840 000000 000000 000000 000000
000000
Moebius: 0
Sum of squares: a^2 + b^2
a = 3514 776401 986872 174070 733209 129673 327241 950873 673372 369609
965291 102998 109899 599898 686750 536018 664732 148375 711432 438199 315006
457855 854921 632037 902485 050909 261792
b = 3514 776401 986872 174070 733209 129673 327241 950873 673372 369609
965291 102998 109899 599898 686750 536018 664732 148375 711432 438199 315006
457855 854921 632037 902485 050909 261792

The bug is not only with this number but I have seen it before. In this particular case it seemts the hiccup is because the number being factored contains a large square.

I belive that YAFU finds a perfect square in the composite being restarted and factors the sqrt of that number, which is fine except if this is the case the number of each prime factor found should be double.

[alexhiggins732]

Also, would it be possible that when the factorization contains multiples of a prime factor the program outputs the exponent instead of the same exponent multiple times?

The problem is say the factorizations contains 2^128, then 2 is printed as factor one 128 seperate lines. There have been cases where so many factors have been displayed I could not scroll up far enough get all the factors.

Yes they are printed to the log file but....

[bsquared]

I see what happens... fermat splits the square into two identical composite factors but only one of the composite factors is re-factored. Trial division removed all of the small squares first, so that's why they were reported correctly.

The easiest fix would be to send both composite factors down the re-factorization path, but that would duplicate a lot of work. I suppose the smart approach would be to test for squaritude first and, if square and composite, just double all of the factors' exponents. I should really have generic ispower test, come to think of it, and use that early in the factor() sequence.

Thanks for the report... and the feature request.

[alexhiggins732]

Originally Posted by bsquared

The easiest fix would be to send both composite factors down the re-factorization path, but that would duplicate a lot of work. I suppose the smart approach would be to test for squaritude first and, if square and composite, just double all of the factors' exponents. I should really have generic ispower test, come to think of it, and use that early in the factor() sequence.

Thanks for the report... and the feature request.

Thank you for your time and effort on this project.

I thought the YAFU did a check for a square already. I remember seeing such a check in the source code.

In my own app, which is nowhere on the level of yours, I check if N is a square then factor the square root after which each exponent is doubled.

Me.Status = ProbablePrime.IsProbablePrime(N)
If Status = PrimeStatus.C Then
Dim root AsBigInteger = Sqrt(N)
If root * root = N Then
Console.WriteLine()
Console.WriteLine(" N is a perfect square. Factoring SQRT(N) = {0} ({1} digits - {2} bits)", root, Digits(root), (BitLength(root)))
log.WriteLine()
log.WriteLine(" N is a perfect square. Factoring SQRT(N) = {0} ({1} digits - {2} bits)", root, Digits(root), (BitLength(root)))
log.WriteLine()
Console.WriteLine()
Me.Add(root, FactoringAlgorithm.U)
ForEach f AsBigIntegerIn Factors.Keys
Factors(f).Exponent *= 2
Next
Else
Me.Add(N, FactoringAlgorithm.U)
EndIf
Else
Me.AddFactor(N, Status, FactoringAlgorithm.P)
EndIf

Where Add() checks compositeness and calls the factoring routine. Now that you mention it, checking for perfect powers should be added.

--------
Edit....

I found your square root check, MPQS.C line 220, v1.20.

[jasonp]

You can add a power detection test, and that will fix the problem here, but may not fix the more general problem: when reducing a number to its prime factorization, you may find a factor that appears in factors that were previously found.

The only way to solve the general version of this problem is to maintain a list of factors found so far, and when a new factor is found you go through the list and divide GCD(new_factor, factor_i) out of composite factor i, then add the GCD result, the smaller old factor and the smaller new factor back into the list. It's easiest to do that recursively.

With the list in hand, you can ask for each composite in the list, then perform a power test and QS or NFS, with new factors added back into the list. Proceed until the list only contains primes. The bookkeeping is easier if you then manually determine how many times each prime in the list divides the input, so there's always only one copy of each prime in the list.

Msieve uses a factor_list_t to do this work; the code using it is in common/driver.c, and it took a surprising amount of fiddling to automatically break down highly composite numbers.

[bsquared]

Originally Posted by alexhiggins732

--------
Edit....

I found your square root check, MPQS.C line 220, v1.20.

Yes, but mpqs was never called for this job... factor() never got that far.

What I need to add is a square/power check in factor().

-----------
Edit: You're welcome! Glad you like it. People that take the time to provide helpful feedback like this is what makes things better, and keeps things fun. :)

[bsquared]

Originally Posted by jasonp

You can add a power detection test, and that will fix the problem here, but may not fix the more general problem: when reducing a number to its prime factorization, you may find a factor that appears in factors that were previously found.

The only way to solve the general version of this problem is to maintain a list of factors found so far, and when a new factor is found you go through the list and divide GCD(new_factor, factor_i) out of composite factor i, then add the GCD result, the smaller old factor and the smaller new factor back into the list. It's easiest to do that recursively.

With the list in hand, you can ask for each composite in the list, then perform a power test and QS or NFS, with new factors added back into the list. Proceed until the list only contains primes. The bookkeeping is easier if you then manually determine how many times each prime in the list divides the input, so there's always only one copy of each prime in the list.

Msieve uses a factor_list_t to do this work; the code using it is in common/driver.c, and it took a surprising amount of fiddling to automatically break down highly composite numbers.

Yep, I've fiddled quite a bit already but it's evident that more is needed ;)

I've got a recursive structure with global factor list already in place - just need to be more thorough about recursing through all composites in the list.

[Mr. Odd]

I upgraded to 1.24 and all of a sudden I get "gnfs has been disabled". Whahappen? It worked fine up until now.

[lorgix]

Same here - 'gnfs has been disabled'

[bsquared]

Sorry, I'd forgotten that NFS wasn't enabled by default. I'll upload a patched 1.24 tonight.

Thanks for the report.