CODEWITHSUNDEEP
algorithmmath
frodo
frodo

Reputation: 1571

Precision issues in finding number of ways to represent a number as a sum of squares

I was trying this problem, http://www.spoj.pl/problems/TWOSQ/ . We have to find the number of different ways to express a number(as large as 10^15) as a sum of squares(without counting twice i.e 5^2 + 1^2 and 1^2 + 5^2 are the same). I have seen this task before and this was how I solved it earlier too. I keep getting Wrong Answer on the judge. Could some one tell me why? or suggest a differnt approach altogether. I have added comments as necessary for understanding . Thanks in advance!.

#include<cstdio>
#include<iostream>
#include<cmath>
using namespace std;

int main()
{
long long X;
cin >> X;
const double EPS = 1e-6;
long long int count = 0;
// add EPS to avoid flooring x.99999 to x
for (int a = 0; a <= sqrt(X/2) + EPS; a++)
{
    long long int b2 = X - a*a; // b^2
    long long int b = (long long int) (sqrt(b2) + EPS);
    if (abs(b - sqrt(b2)) < EPS) // check b is an integer
        count++;
}
cout << count << endl;

}

Upvotes: 1

Views: 210

Answers (4)

frodo
frodo

Reputation: 1571

I tried quite a bit but was still getting a Wrong Answer for one reason or another. Hence I looked at a different approach and was able to figure this pretty fast solution out.

 long long L=res=0, R=(long long)sqrt(X) + 1;     
while (L<=R)
{
    long long Y = L*L + R*R;
    if (Y > X) R--;
    else if (Y < X) L++;
    else res++, L++;
}

Upvotes: 0

Samy Vilar
Samy Vilar

Reputation: 11130

I think I know what you are trying to do, its rather cool ... try this:

#include<cstdio>
#include<iostream>
#include<cmath>

using namespace std;

typedef unsigned long long int data_type;

int main()
{
    data_type value = 0,
              count = 0, 
              index = 0, 
              max = 0,         
              b_1 = 0; 

    cin >> value;
    max = (data_type)sqrt(value);
    char *flags = new char[max];
    memset(flags, 0, max);

    for (index = 0; index < max; index++)
    {   
        b_1 = value - (index * index);        

        if (!((data_type)sqrt(b_1) - sqrt(b_1)) && !flags[(data_type)(sqrt(b_1))] && !flags[index])
        {
            flags[(data_type)sqrt(b_1)] = 1;
            flags[index] = 1;
            count++;                          
        }
    }        

    cout << count << endl;
}

as far as I can tell you square then take the difference, the check if that difference itself is a perfect square by taking the square root and checking for any residual floating point values, kind of cool.

This hasn't being thoroughly tested As such if theres any bugs or issues, I apologize.

As I thought its nice trick the problem is repetitions.

Note that (10**15)**.5 is roughly 32 megs so thats how much it will allocate
good luck.

Upvotes: 0

Jeffrey Sax
Jeffrey Sax

Reputation: 10313

There are two problems that I can see.

  1. In your calculation of b2, you use the expression a*a. If a is just an int, then this will overflow very quickly.
  2. Your value for EPS is way too large. You will get false positives.

Double precision floating-point numbers are stored using up to 53 significant bits. This means that all integers up to about 8e15 can be represented exactly. For the square root of such a number to be rounded correctly, you'd need about double the precision, so that leaves you with 4e15, still within your range.

So, I would do two things:

  1. Change all my variables to doubles.
  2. Do away with EPS entirely and use exact comparisons. They should work fine within the range you specify (up to X = 1e15).

Upvotes: 1

Erick Wong
Erick Wong

Reputation: 165

Your use of EPS in the main loop is fundamentally flawed (it is conceivable that it could escape any actual errors for Diophantine reasons, but this would require more thought than is warranted). Say X/2 is very close to a perfect square, for instance 10^14-1. Then sqrt(X/2) will be about 10^7 - 0.5*10^-7, which will underflow a double precision float and round to 10^7.

Why not just stop the loop when 2*a*a > X? And of course, fix your integer overflow hazards.

Upvotes: 0

Related Questions