Variable might be set by a nonscalar (three variables in pythagorean triplet)

Hi everyone,

Im working for project euler and I try to do #9 (Special Pythagorean Triplet). I have encountered couples of errors, one of which is something that occurs from my trying to make 3 arrays probably, some of them are named Variable might be set by a nonscalar and other errors which I don't understand.

function y = pythtrip(a,b,c)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here
y = a*b*c;

Unrecognized function or variable 'a'.

a = 0:750;
b = 0:750;
c = 0:750;
while (a>b) && b>c && (a^2)+(b^2)==(c^2) && a+b+c == 1000;
    pythtrip(a,b,c)
y;
end

[Torsten]

Before starting in your project, you should visit

https://matlabacademy.mathworks.com/details/matlab-onramp/gettingstarted

The course is free of costs and helps you to learn the basics of the new software.

[Konstantin]

Though, I know the basics, except things that i had forgotten.

Is it difficult to tell me couple of things to improve ?

[Torsten]

I guess you want to find all Pythagorean triples (a,b,c) with a < b < c and a + b + c = 1000.

ifound = 0;
for a = 1:333         %  a < b < c and a + b + c = 1000 --> a <= 333
    lb = a + 1;
    ub = 500 - a/2 - 1 + 0.5*mod(a,2);
    for b = lb:ub   % a < b and c = 1000 - a - b > b --> a < b & b < 500 - a/2
        c = 1000 - a - b;
        if a^2 + b^2 == c^2
            ifound = ifound + 1;
            A(ifound) = a;
            B(ifound) = b;
            C(ifound) = c;
        end
    end
end
ifound

ifound = 1

A

A = 200

B

B = 375

C

C = 425

[Konstantin]

Yes, that's it. ( https://projecteuler.net/problem=9 ) Thank you

It's the smart way to define array for only one variable and define other in relationship with the first one.

[Torsten]

Optimized search code according to

https://stackoverflow.com/questions/61755082/pythagorean-triplet-with-given-sum

s = 1234571816;
ifound = 0;
% Solve (1) b + c = s - a, (2) c^2 - b^2 = a^2 for b as a function of a
% Result : b = s*(s-2*a)/(2*(s-a)) 
% Consequence: Only a loop over a is required, b and c follow 
for a = 1:s/3
  num = s*(s-2*a);   
  den = 2*(s-a);      
  if mod(num,den)==0   
    b = num/den;
    c = s - a - b;
    if a < b
      ifound = ifound + 1;
      A(ifound) = a;
      B(ifound) = b;
      C(ifound) = c;
    end
  end
end
ifound

ifound = 5

A

A = 1×5

8 72621864 116469032 258255456 268152008

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

B

B = 1×5

617285904 578705543 552985297 454001280 446007904

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

C

C = 1×5

617285904 583244409 565117487 522315080 520411904

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[John D'Errico]

Note: even the code provided by @Torsten will fail for only a little bit larger values of the sum. In the example he shows, the loop will run up to

1234571816/3

ans = 4.1152e+08

So 400 million iterations through that for loop. The solution I provided requires only the ability to factor the sum.

s1 = 1234571816;
timeit(@() factor(s1))

ans = 2.5282e-04

Now try the same thing, but for a larger sum.

s2 = sym('1234634647253636364545');
factor(s2)

ans = 

timeit(@() factor(s2))

ans = 0.0081

However, you might want to catch up on your reading while you wait for a loop from 1 to s2/3 to terminate.

[Torsten]

I think OP would have been happy already with his three-loops solution: a good horse does not jump higher than it has to :-) It's sad somehow.

[John D'Errico]

@Torsten

But you are making EXACTLY my point. If you are solving PE problems, then a triply looped solution will stop being viable very quickly. Even an O(S) solution, as you showed, is useless if S becomes just a little larger, say 1e20. And PE problems are notorious for exactly that, at least once you get past the first dozen or so. Then you need to be looking for better solutions, not brute force. And even a linear solution can be poor at some level.

[Torsten]

Yes, I agree: your solution is better than mine (or since it's not mine: the one that I found as classified "optimal" because of complexity O(S)). Maybe it's more difficult to put your method into formal MATLAB code where you only have to specify the value of the sum and all the triples are automatically generated. But in principle it's more efficient.

Maybe your method is hard to find since the complexity of factoring is not known to be polynomial in the input size.

Answer 1

[John D'Errico]

Let me suggest that brute force is NEVER the good way to solve a project Euler problem. NEVER. When brute force is applied, most of brute froce codes you would write would still be running long after both you and your computer are piles of dust. (Ok, I'll admit that the Sudoku PE problem was doable by me merely solving each sudoku puzzle by hand, then computing a result. That was ok, since I like doing sudoku puzzles. But that was one of the very, very few PE problems I solved on paper, and I recall having solved around 300 of them almost entirely in MATLAB.) I will comcede that this problem is probably doable via brute force, as a pretty low numbered problem. But that still makes it a poor choice of method.

Instead, start with the Euclidean formula for Pythagorean triples. ALWAYS do some reading and thinking BEFORE you try to write code for a PE problem.

https://en.wikipedia.org/wiki/Pythagorean_triple

In there, you will learn how to represent all Pythagorean triples as the triple:

{a,b,c} = {k*(m^2-n^2), 2*k*m*n, k*(m^2+n*2)}

In that formula, we know that m>n is necessary. What do you know about the desired triples? Their sum must be 1000, What can you get from that formula?

a+b+c == 2*k*m*2 + 2*k*m*n == 2*k*m*(m+n) == 1000

Now look for all such values of k,m,n, where m>n, such that the produt holds. It looks like solution for a sum as small as 1000 can be done using pencil and paper. And since they tell you that only one solution exists, you can do it via some pretty fast trial and error, just by looking at the factor pairs of a number.

factor(1000)

ans = 1×6

2 2 2 5 5 5

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

We know that one of those factors of 2 is already taken.

factorpairs = @(N) [divisors(N)',N./divisors(N)'];

First, suppose k is 1, then we can find M and N

fp = factorpairs(500)

fp = 12×2

1 500 2 250 4 125 5 100 10 50 20 25 25 20 50 10 100 5 125 4 250 2 500 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Each such pair of factors would imply a value for m and n. Remembering that m>n is a requirement for the Euclid formula to work, we might do this:

mnsolve = @(fp) [fp(:,1), fp(:,2) - fp(:,1)];
mnsolve(fp)

ans = 12×2

1 499 2 248 4 121 5 95 10 40 20 5 25 -5 50 -40 100 -95 125 -121 250 -248 500 -499

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

And when k==1, only one such pair seems to yield a viable pythagorean triple.from that set. (Remember the rules about m and n. They must be both positive, and must have m>n.) You might decide to test it to see if that pair indeed satisfies your needs.

triplefun = @(k,m,n) k.*[m.^2-n.^2, 2*m.*n, m.^2+n.^2)];

Now you could try it again for k=2, which would require the factor pairs of 250.

Honestly, I don't recall my solution method for this problem, but I would guess I did something fairly simple like this. I could look it up, since I think I saved the MATLAB codes for all of those PE problems.

[John D'Errico]

Just for kicks, using the ideas in my answer, how would I find a pythagorean triple with a designated sum a+b+c=123456?

S = 123456;
factor(S)

ans = 1×8

2 2 2 2 2 2 3 643

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

factorpairs = @(N) [divisors(N)',N./divisors(N)'];
mnsolve = @(fp) [fp(:,1), fp(:,2) - fp(:,1)];

% with k == 1:
k = 1;
fp = factorpairs(S/2/k)

fp = 24×2

1 61728 2 30864 3 20576 4 15432 6 10288 8 7716 12 5144 16 3858 24 2572 32 1929 48 1286 96 643 643 96 1286 48 1929 32

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

mnsolve(fp)

ans = 24×2

1 61727 2 30862 3 20573 4 15428 6 10282 8 7708 12 5132 16 3842 24 2548 32 1897 48 1238 96 547 643 -547 1286 -1238 1929 -1897

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

There are clearly no valid solutions when k==1. But, when k==643 (there may be other solutions, I did not look further through all of the divisors of S to learn if there may be others. I'd not be amazed if there were.)

k = 643;
fp = factorpairs(S/2/k)

fp = 12×2

1 96 2 48 3 32 4 24 6 16 8 12 12 8 16 6 24 4 32 3 48 2 96 1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

mnsolve(fp)

ans = 12×2

1 95 2 46 3 29 4 20 6 10 8 4 12 -4 16 -10 24 -20 32 -29 48 -46 96 -95

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

We do see there parameters m and n that will generate a triple. m>n is the requirement, and both must be positive. I could have automated this fairly easily, but the explanations are I think more useful.

triplefun = @(k,m,n) k.*[m.^2-n.^2, 2*m.*n, m.^2+n.^2];
T = triplefun(643,8,4)

T = 1×3

30864 41152 51440

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

sum(T)

ans = 123456

Is this a pythagorean triple?

dot(T.^2,[1 1 -1])

ans = 0

As expected, indeed it is. Note that the triple I did find is a pretty large one. And brute force loops would have taken a fair amount of time to produce a result. As well, you should see the above computation took essentially no CPU time to find. What if the sum had been given as something on the order of 1e12 instead? What if it had been as large as 1e100? That would be a problem more worthy of Project Euler. As long as I used syms in the computations there to avoid double precision issues, it would still have taken virtually no effort to solve. Again, this is what you need to take away from Project Euler, to look more deeply into a problem than just brute force. While brute force can work for small problems, it often becomes intractable, and that can happen very quickly.