CODEWITHSUNDEEP
c++
E.O.
E.O.

Reputation: 814

Calculating pythagorean triples

I was writing some code to calculate Pythagorean triples, but then I got some values which were not solutions. Here is the code:

#include <iostream>
#include <math.h>

using namespace std;

int main()
{
int j, counter = 0;
double candidate, sqrnbr;
for (int i = 400; i <= 500; i++) //note that the range is only from 400-500 since I am still testing the code
{
    for (j = i; j <= 500; j++)
    {
        candidate = i*i+j*j;
        sqrnbr=sqrt(candidate);
        if (candidate/sqrnbr==sqrnbr)
        {
           cout << i << " and " << j << " and " << sqrnbr << endl;
           counter++;
        }
    }

}
cout << counter << endl;
system("pause");
return 0;
}

One of the outputs it gave was 408, 410 and 578.415. Does anyone know why this happens?

Upvotes: 2

Views: 4837

Answers (6)

Paul Baxter
Paul Baxter

Reputation: 1175

This program generates primitive triplets. It takes a command line argument for the start and end. It also uses 8 worker threads for large numbers. I have tested this with generating primitives so a b and c < 5,000,000. 

#include <cmath>
#include <iostream>
#include <thread>
#include <fstream>
#include <chrono>

using namespace std;

#define MAX_INT64       (0X7FFFFFFFFFFFFFFFLL)
#define MOD(a,b)        (a%b)
#define SQRT(x)         (sqrt(x))
#define TRUNC(x)        ((MYTYPE)(x))
#define NUM_BLOCKS      (1000)
#define WORKER_THREADS  (8)

typedef long long MYTYPE;
typedef long double MYFLOATTYPE;

// Structure for a Triple
typedef struct Triple
{
    MYTYPE A;
    MYTYPE B;
    MYTYPE C;
} Triple;

// Structure for a block of Triples
typedef struct TripleBlock
{
    Triple* Triples;
    long NumTriples;
    long NumTriplesAllocated;
    TripleBlock* Next;
} TripleBlock;

// Structure for worker thread
typedef struct ThreadStartData
{
    MYTYPE Min;
    MYTYPE StartMax;
    MYTYPE Max;
    MYTYPE NumDone;
    MYTYPE NumThreads;
    MYTYPE Complete;
    TripleBlock* Block;
} ThreadStartData;

void PrintTriples(ThreadStartData* lpParam, ofstream* file);
TripleBlock* AllocateTripleBlock(long numtriples);
void AddTriple(MYTYPE a, MYTYPE b, MYTYPE c, TripleBlock* block);
void CalculateTriples(ThreadStartData* lpParam);
void ProgressTriples(ThreadStartData* lpParam);

int main(int argc, char* argv[])
{
    MYTYPE min = 2;
    MYTYPE max = 0;

    if (argc == 1)
    {
        max = MAX_INT64;
    }
    else if (argc == 2)
    {
        max = atol(argv[1]);
    }
    else if (argc == 3)
    {
        min = atol(argv[1]);
        max = atol(argv[2]);
    }
    else
    {
        cout << argv[0] << " [min] max" << endl;
        return 0;
    }
    if (max <= min || min < 1)
    {
        cout << "Invalid parameters" << endl;
        return 0;
    }

    int numThreads = WORKER_THREADS;
    if ((max - min) <= (WORKER_THREADS * 2))
        numThreads = 1;

    MYFLOATTYPE inc = ((MYFLOATTYPE)max - (MYFLOATTYPE)min) / ((MYFLOATTYPE)numThreads);
    MYFLOATTYPE current = (MYFLOATTYPE)min;

    // set up ThreadStartData
    ThreadStartData* data = new ThreadStartData[numThreads];
    for (int i = 0; i < numThreads; i++)
    {
        data[i].Min = min;
        data[i].Max = max;
        data[i].NumDone = 0;
        data[i].NumThreads = numThreads;
        data[i].Complete = 0;
        data[i].Block = AllocateTripleBlock(NUM_BLOCKS);

        current += inc;
        MYTYPE startmax = TRUNC(current);
        if (startmax > max)
            startmax = max;
        if (i == numThreads - 1)
        {
            startmax = max;
        }

        data[i].StartMax = startmax;

        if (data[i].StartMax > max)
            data[i].StartMax = max;

        min = startmax + 1;
        if (min > max)
            min = max;
    }

    // start the threads
    thread* hThreadArray = new thread[numThreads];
    for (int i = 0; i < numThreads; i++)
    {
        hThreadArray[i] = thread(CalculateTriples, &data[i]);
    }
    thread hProgressthread = thread(ProgressTriples, data);

    // wait for threads to complete
    for (int i = 0; i < numThreads; i++)
    {
        hThreadArray[i].join();
    }
    hProgressthread.join();

    // print the results
    ofstream* myfile = new ofstream;
    myfile->open("triples.txt");
    PrintTriples(data, myfile);
    myfile->close();

    // cleanup
    for (int i = 0; i < numThreads; i++)
    {
        TripleBlock* block = data[i].Block;
        while (block)
        {
            TripleBlock* next = block->Next;
            if (block->NumTriplesAllocated > 0)
                delete[] block->Triples;
            delete block;
            block = next;
        }
    }

    delete[] data;
    delete[] hThreadArray;

    return 0;
}

/// <summary>
/// Progresses thread for the triples.
/// </summary>
/// <param name="data">The data.</param>
void ProgressTriples(ThreadStartData* data)
{
    MYTYPE numThreads = data->NumThreads;
    MYTYPE max = data->Max;
    MYTYPE lastPercent = -1;
    bool done = true;
    ThreadStartData* dataStart = data;

    do
    {
        MYTYPE total = 0;
        done = true;
        data = dataStart;

        /// total up the number done
        for (int i = 0; i < numThreads; i++, data++)
        {
            total += data->NumDone;

            // check if thread is done
            if (data->Complete == 0)
                done = false;
        }

        // see if all threads are completed
        if (done)
        {
            cout << "done." << endl;
            break;
        }

        MYFLOATTYPE percentdone = ((MYFLOATTYPE)total / (MYFLOATTYPE)max) * 100;
        if (lastPercent != TRUNC(percentdone))
        {
            lastPercent = TRUNC(percentdone);
            cout << lastPercent << "% done.  (" << total << " of " << max << ")" << endl;

            this_thread::sleep_for(chrono::seconds(2));
        }
    } while (!done);
}

/// <summary>
/// Calculates the triples.
/// </summary>
/// <param name="data">The data.</param>
void CalculateTriples(ThreadStartData* data)
{
    MYTYPE min = data->Min;
    MYTYPE max = data->Max;
    MYTYPE startmax = data->StartMax;
    TripleBlock* block = data->Block;

    for (MYTYPE a = min; a >= min && a <= startmax; a++)
    {
        data->NumDone++;

        MYFLOATTYPE asquared = (MYFLOATTYPE)a * (MYFLOATTYPE)a;
        for (MYTYPE b = a + 1; b > a; b += 2)
        {
            MYFLOATTYPE bSquared = (MYFLOATTYPE)b * (MYFLOATTYPE)b;
            MYFLOATTYPE cSquared = asquared + bSquared;
            MYFLOATTYPE cSquareRoot = SQRT(cSquared);
            MYTYPE c = TRUNC(cSquareRoot);

            // check for max and overflow
            if (c > max || c < b)
                break;

            // check for a true square
            if (cSquareRoot != c)
                continue;

            // check gcd
            MYTYPE gcd = a;
            for (; gcd >= 1; gcd--)
            {
                if ((MOD(c, gcd) != 0) || (MOD(b, gcd) != 0) || (MOD(a, gcd) != 0))
                    continue;

                break;
            }

            // if gcd != 1 then no triplet
            if (gcd != 1)
            {
                continue;
            }

            AddTriple(a, b, c, block);
            break;
        }
    }

    data->Complete = -1;
}

/// <summary>
/// Adds the triple.
/// </summary>
/// <param name="a">a.</param>
/// <param name="b">b.</param>
/// <param name="c">c.</param>
/// <param name="block">block.</param>
void AddTriple(MYTYPE a, MYTYPE b, MYTYPE c, TripleBlock* block)
{
    while (block->NumTriples >= block->NumTriplesAllocated)
    {
        if (block->Next == NULL)
        {
            block->Next = AllocateTripleBlock(NUM_BLOCKS);
        }
        block = block->Next;
    }

    Triple* triple = &(block->Triples[block->NumTriples]);
    triple->A = a;
    triple->B = b;
    triple->C = c;
    block->NumTriples++;
}

/// <summary>
/// Allocates a triple block.
/// </summary>
/// <param name="numtriples">The numtriples.</param>
/// <returns>TripleBlock *</returns>
TripleBlock* AllocateTripleBlock(long numtriples)
{
    TripleBlock* block = new TripleBlock;
    memset(block, 0, sizeof(TripleBlock));

    block->Triples = new Triple[numtriples];
    memset(block->Triples, 0, numtriples * sizeof(Triple));
    block->NumTriplesAllocated = numtriples;

    return block;
}

/// <summary>
/// Prints the triples.
/// </summary>
/// <param name="data">The data.</param>
/// <param name="myfile">The file.</param>
void PrintTriples(ThreadStartData* data, ofstream* theFile)
{
    for (int i = 0; i < data->NumThreads; i++)
    {
        TripleBlock* block = data[i].Block;
        while (block)
        {
            Triple* triple = block->Triples;
            for (MYTYPE i = 0; i < block->NumTriples; i++, triple++)
            {
                (*theFile) << "  " << triple->A << ", " << triple->B << ", " << triple->C << endl;
            }
            block = block->Next;
        }
    }
}

Upvotes: 1

reza
reza

Reputation: 1

#include <iostream>

using namespace std;

int main()
{
   int side1, side2, side3;
   int counter=0;
   for(side1=1;side1<=500;side1++)
    for(side2=1;side2<=500;side2++)
     for(side3=1;side3<=500;side3++)
       {
        int t1, t2, t3;
        t3=side1 * side1 + side2 *side2;
        t2=side1 * side1 + side3 * side3;
        t1=side2 * side2 + side3 * side3;
        if(t3==side3 * side3 || t2==side2 * side2 || t1==side1 * side1)
        {
            cout<<side1<<" , "<<side2<<", "<<side3<<endl;
            counter++;

        }
       }
   cout<<endl<<counter/6;
}

Upvotes: 0

Nic Foster
Nic Foster

Reputation: 2894

Your algorithm is incorrect.

Try changing your if to:

if ( (candidate / sqrnbr) == floor(sqrnbr))

This should make sure you end up inside your if statement when you have whole numbers.

Upvotes: 0

Sergey Kalinichenko
Sergey Kalinichenko

Reputation: 726479

Comparing doubles for equality is a bad idea in general. You should make candidate an int, and change the program around like this:

#include <iostream>
#include <math.h>

using namespace std;

int main() {
    int j, k, counter = 0, candidate;
    double sqrnbr;
    for (int i = 400; i <= 500; i++) {
        for (j = i; j <= 500; j++) {
            candidate = i*i+j*j;
            sqrnbr=floor(sqrt(candidate));
            if ((sqrnbr*sqrnbr)==candidate) {
               cout << i << " and " << j << " and " << sqrnbr << endl;
               counter++;
            }
        }
    }
    cout << counter << endl;
    system("pause");
    return 0;
}

Note the floor around the sqrt: it drops the non-integer part of the square root.

Upvotes: 2

Blender
Blender

Reputation: 298046

Your logic was a bit flawed.

I reformatted your code and it works for me:

#include <iostream>
#include <math.h>

using namespace std;

int main() {
  int counter = 0;
  float candidate;

  for (int i = 1; i <= 500; i++) {
    for (int j = 1; j <= 500; j++) {
      candidate = sqrt(i*i + j*j);

      if (int(candidate) == candidate) {
        cout << i << " and " << j << " and " << candidate << endl;
        counter++;
      }
    }
  }

  cout << counter << endl;
  return 0;
}

To check whether sqrt(sum) is an int, just check whether casting it to an int changes its value.

Upvotes: 1

Daniel Fischer
Daniel Fischer

Reputation: 183873

You calculate a double and never check whether it's an integer value, so you get double values. You don't get a lot more because the square root is generally not representable as a double, so

candidate/sqrnbr == sqrnbr

is false for many values of candidate.

To generate Pythagorean triplets, use only integer arithmetic. And there are better methods than brute force, there's a classic formula, all Pythagorean triplets are of the form

d*(m^2 - n^2), d*2*m*n, d*(m^2 + n^2)

for some d, m, n.

Upvotes: 4

Related Questions