[GAP Forum] Different TwoSquares algorithm for large input?

[Jason B Hill]

The available documentation helps with the original question a bit. I'll
mention some options in a bit. First, notice the difference between gap and
sage.

Sage's 'two_squares()' method indeed calls gap's TwoSquares function, and is
literally written as "gap.gap.eval('TwoSquares(%s)'%n)". The documentation
claims that two_squares(n) will write n as the sum of two squares if
possible, and return a value error otherwise.

It should be noted that gap's TwoSquares is more complex, as found in 15.6
of the gap reference manual: "TwoSquares returns a list of two integers x\le
y such that the sum of the squares of x and y is equal to the nonnegative
integer n... If no such representation exists Two Squares will return fail.
TwoSquares will return a representation for which the gcd of x and y is as
small as possible. It is not specified which representation TwoSquares
returns if there is more than one."

The sage method was clearly written with the idea that it may expanded by
continued fractions (via pari) in the future. But, for right now, the sage
documentation isn't exactly precise and I'll submit a patch to clear this
up. The two main concerns are (a) that the value exception raised looks more
like an error (someone may wish to test if n has a two squares
representation and if not then this method doesn't make that entirely clear)
and (b) the documentation does not explain that a lowest gcd representation
is returned. It would also be nice to have (1) a potentially faster method
that returns any two squares representation regardless of gcd considerations
and (2) a boolean method that states if a two squares representation exists
(without calculating such).


Some options for making your code faster/better:

1) Factoring large integers in sage will be slower than in gp/pari (even
though 'factor()' in sage does use gp/pari), as the gp-prime-limit is fixed
and sage forces a global proof attribute on factorizations. Since you have
sage installed, you might as well use {{{sage -gp -p prime_limit}}} with a
high prime limit (and maybe more memory?) to factor anything too serious.

2) To test if a two squares representation exists, divide from n all of the
primes in the factorization of n not of the form 4k+3. The remaining product
is a square iff n has a two squares representation. This is probably how gap
determines if it should return "fail".

3) If you call a gap or gp command from sage that will return something like
a list or list of lists, it is better to cast it back as a sage object:
sage: gap('FactorsTD(12^25+25^12)').sage()
Objects like the following will obviously require more parsing:
sage: gp('factor(12^25+25^12)')
[13, 1; 19, 1; 727, 1; 35149, 1; 151142573749569397, 1]



Jason