Inefficient square root

[GregoryMorse]

Description
During the square root in GNFSCore/SquareRoot/SquareFinder.cs, it will find the positive and negative roots via inverse modulo. Then it compares it to the positive and negative (inverse) original polynomial which does not actually do anything nor make sense. Finally it does not make use of this anyway, and uses a cartesian product and searches all 2^d values.

According to the thesis An Introduction to the General Number Field Sieve by Matthew E. Briggs cited here, the norm function should be used to disambiguate this. It is computed as alpha^((p^d-1)/(p-1) where alpha should be both of the potential roots. However I also cannot figure out how to get this to work, as it is unclear what is a positive versus negative norm. The key is being consistent across the splitting field. But trial and error is an oversight that should be correctable. Unfortunately I cannot find the details to get this working.

[AdamWhiteHat]

tl; dr: Yeah, you are correct; Something is broken with the final square root step. I am looking at that now. I spent several hours on it last night, and I plan to spend all of today on it as well. Ill probably keep plugging away at it until I get it working. There is some progress being made. Ive identified a problem area, due to your comments. I'm currently figuring out what to do about it and also busy looking through the literature and verifying my implementation in certain areas.

Thanks for bringing this to my attention. I was aware of it, I was just demotivated to find the issue, but you submitting a bug report and me knowing that other people are looking at it has motivated me once again. Well, and I have been meaning to take a look at this again. So thanks for the push.

It was working at one point in time, and then * something * changed and now the difference of two squares fails to yield non-trivial factors.
The code was under source control at the time, but I have scoured the commits but its not clear what changed to cause the issue. Any changes to the square root code was heavily scrutinized and I could not see the problem, so it may be that a change that lies elsewhere is the culprit, such as in the relations factorization or the matrix step.

You are correct that taking the cartesian product and trying every possible combination is wrong. This is an artifact left over from when I was flailing, trying to make the square root step produce factors again--evidence of my desperation at the time. I probably left it in because it works sometimes to find factors, although if it does its not by any principled approach or theory, its either brute forcing it or it could very well be that it does so entirely due to a numerical coincidence, and the fact that the reference/test numbers being factored are so small.

I'm not too sure what you mean by "positive and negative roots via inverse modulo", but by the fact that you say inverse modulo and say that I'm comparing them, I take this to mean in a method called AlgebraicSquareRoot, how I compare the startPolynomial with the resultSquared1 and the modulo inverse of the start polynomial with resultSquared1, looking for a match?

Yeah, you know, looking at it, that is odd and I'm not sure that's right. The first condition should be enough, where the square root polynomial found by the Tonelli-Shanks algorithm method (which I'm not calling Tonelli-Shanks , I'm calling FiniteFieldArithmetic.SquareRoot for some reason) and the modular inverse of that polynomial, when squared, should both be the same.

At any rate, I'm looking into it. I was able to get it working using Matthew E. Briggs' paper once, I'm sure I can get it working again.
Ill keep you informed.

[GregoryMorse]

So that thesis seems to use the idea here: https://www.math.u-bordeaux.fr/~jcouveig/publi/Cou94-2.pdf

Taking the inverse of startPolynomial doesn't make any sense. It will never match.

Basically you can compute N(alpha) and Np(alpha)=alpha^((p^d-1)/(p-1)). N(alpha) isn't computed from your potentially wrongly signed square root.

The Briggs paper handwaves that important detail.

It's N(alpha)=N(f'(alpha))p1^fp1p2^fp2

Basically the result from the sieve as a square should be decomposed. In the same way as the exponents were determined when building the X matrix. Ord_p(N(a+b*alpha)). They are guaranteed to be even powers as the result is a square.

At least this is my interpretation. Also working on it. Interesting to hear your thoughts about it. Good examples are scarce for some of these detailed steps.

[GregoryMorse]

Yea that solves it...

So when you are checking that prod(p)==(-b)^d*f(-b/a) you should save all the first-order p values!

Why? Because you can tabulate them after the fact, all of them have degree 1 by definition there, but if you count each prime, they will all end up with even powers after getting the result vector - you dont just want to compute prod(a+b*alpha)=beta^2 but also the order of each prime. For example I made a dictionary like so:

{2: 2, 67: 4, 7: 2, 31: 2, 61: 2, 73: 4, 41: 2, 89: 4, 79: 2, 23: 2, 43: 2, 29: 2}

Since my smooth pairs (a,b) found as the result were: { (-2, 1),(1, 1),(13, 1),(15, 1),(23, 1),(3, 2),(33, 2),(5, 3), (19, 4),(14, 5),(15, 7),(119, 11),(175, 13),(-1, 19),(49, 26) }

And corresponding to those smooth pairs were the following primes found during the smooth pair finding process:
[{2, 67}, {7}, {31}, {61, 7}, {73, 67}, {41}, {89, 79}, {89}, {89, 73}, {2, 61, 23}, {73, 23}, {89, 41, 67}, {43, 67, 79}, {73, 29, 31}, {43, 29}]

Now all those by definition have to be even, as they are squares. But we can halve them, and use them to compute the norm. Not sure if I explained this well, but this is the technique outlined in the paper. N(f'(alpha))*prod(all the primes to half the power) == root^((p^d-1)/(p-1)) gives the correct root. However, I cannot figure out N(f'(alpha)), this norm is ill defined in the papers. Too much hand waiving, the authors were not thorough.

The limitation to this approach, is the polynomial degree must be odd. Otherwise for even degree, perhaps the Cartesian product is necessary. Of course you dont need a full Cartesian product. You can probably fix one value as all negative or all positive will should both work. Its about consistency. So pivoting the first value would arbitrarily half the search space at least.

[GregoryMorse]

Finally a paper that went through the details: https://infoscience.epfl.ch/record/164498/files/164498.PDF

N(f'(alpha)) can be computed via f'(alpha)^((p^d-1)/(p-1))

Will be interesting to see how you handle it. Not sure if you assume odd degree polynomials, but they can be handled quite efficiently compared to the even counterpart.

[AdamWhiteHat]

The main issue now is lies in step 5.9 (page 58) of Matthew E Briggs paper: Determining Applicable Finite Fields

The issue is the method Polynomial.Field.GCD is defective.
With this not working, the algorithm is unable to tell if f(x) is irreducible over Z/pZ, for some random 'inert' prime p.

IIRC, I really struggled to get this method implemented, and apparently its still not right. Nevermind, I was thinking of the multivariate polynomial GCD, which requires having to calculate the Grobner basis. A single variable polynomial GCD is pretty straight forward. Hopefully it there is just a silly bug I overlooked.

Once this is working correctly, I believe the GNFS algorithm will work, because the algorithm finds the factors for 45113 if it is handed the 3 primes in the worked example that makes applicable finite fields. as demonstrated in the SmallFactorizationTest_45113 test class in the test project.