# HP Forums

Full Version: A Sunday Programming Challenge
You're currently viewing a stripped down version of our content. View the full version with proper formatting.

Project Euler problem 370

Quote:
Let us define a geometric triangle as an integer sided triangle with sides a <= b <= c so that its sides form a geometric progression, i.e. b^2 = a*c.

An example of such a geometric triangle is the triangle with sides a = 144, b = 156 and c = 169.

There are 42 geometric triangles with perimeter <= 100. Appropriately for this problem F(100)=42.

Write a program (42S or platform of your choice) optimized for speed then for size that can solve F(1000) on a real 42S or F(1e8) on PC emulator in a reasonable time.

Edited: 11 Mar 2012, 1:05 a.m.

Unoptimized and very slow, I fear:

```00 { 71-Byte Prgm }
01>LBL 00
02 STO 00
03 3
04 BASE÷
05 1E-3
06 +
07 STO 01
08 IP
09 STO 02
10>LBL 01
11 RCL 01
12 1
13 BASE-
14 STO 03
15>LBL 02
16 RCL 01
17 IP
18 X^2
19 RCL÷ 03
20 STO 04
21 IP
22 RCL 03
23 RCL+ 01
24 X<Y?
25 GTO 04
26 RCL 04
27 FP
28 X!=0?
29 GTO 03
30 RCL 00
31 RCL 04
32 RCL 01
33 BASE+
34 RCL+ 03
35 X>Y?
36 GTO 03
37 1
38 STO+ 02
39>LBL 03
40 1
41 STO- 03
42 GTO 02
43>LBL 04
44 DSE 01
45 GTO 01
46 RCL 02
47 .END.
100 XEQ 00  -->   42     (1 min 19 sec, real HP-42S)
1000 XEQ 00  -->  532     (~1 second (Free42 Decimal)
```

Quote:
There are 42 geometric triangles with perimeter <= 100.

I get a significant different number of 'goemetric integer traiangle' !

I may have miss samething, or badly interpret the conditions.

I count 83 integer sided triangles (a,b,c) with a <= b <=c and perimeter P less to 100 who satisfed b2 = a.c.

Are you sure you are not counting integer sided triangles with different side's length restrictions (such as 1 < a < b < c ) ?

Edited: 11 Mar 2012, 5:45 p.m. after one or more responses were posted

Quote:
I count 83 integer sided triangles (a,b,c) with a <= b <=c and perimeter P less to 100.

Remember it must also satisfy b^2 = a*c

Are you checking that condition also?

What is step "24 X<Y?" testing for? It looks like a+b < c as a loop termination condition. Why is this necessary?

Apart from the trivial cases of a=b=c this program would seem to produce (I converted it to C to play with so could be wrong here) these nine extras:

``` 25  30  36
20  30  45
18  24  32
16  24  36
16  20  25
12  18  27
9  12  16
8  12  18
4   6   9
```

My next question is where are triples like 1 3 9 or 14 28 56? These look to be in geometric progression and satisfy a c = b2. In fact I get fifty of these which include the nine above:

``` 25  30  36
20  30  45
18  30  50
16  28  49
14  28  56
13  26  52
18  24  32
16  24  36
12  24  48
9  24  64
11  22  44
9  21  49
7  21  63
16  20  25
10  20  40
8  20  50
12  18  27
9  18  36
6  18  54
8  16  32
4  16  64
9  15  25
5  15  45
3  15  75
7  14  28
4  14  49
9  12  16
8  12  18
6  12  24
4  12  36
3  12  48
2  12  72
5  10  20
4  10  25
2  10  50
3   9  27
1   9  81
4   8  16
2   8  32
1   8  64
1   7  49
4   6   9
3   6  12
2   6  18
1   6  36
1   5  25
2   4   8
1   4  16
1   3   9
1   2   4
```

Add to this the 33 a=b=c cases and I'm at 83.

What am I missing or is the problem stated incorrectly?

- Pauli

I'm getting 83 too (see above). They all satisfy the condition as stated.

- Pauli

I was getting 83 also, but then I realized we are looking for triangles. In a triangle the length of one side cannot exceed the sum of the lengths of the other two sides, as we know. That's what step 24 is testing for.

Gerson.

Quote:
What is step "24 X<Y?" testing for? It looks like a+b < c as a

How do you know this is a triangle:

```  1   9  81
```

Yes, all the 83 I count satisfy b2=ac condition.

I also get 50 by considering only a < b < c.

```a	b	c	P=a+b+c			a	b	c	P=a+b+c
1	2	4	7	nat		1	1	1	3
1	3	9	13	nat		2	2	2	6
2	4	8	14	nat		3	3	3	9
4	6	9	19			4	4	4	12
1	4	16	21	nat		5	5	5	15
3	6	12	21	nat		6	6	6	18
2	6	18	26	nat		7	7	7	21
4	8	16	28	nat		8	8	8	24
1	5	25	31	nat		9	9	9	27
5	10	20	35	nat		10	10	10	30
9	12	16	37			11	11	11	33
8	12	18	38			12	12	12	36
3	9	27	39	nat		13	13	13	39
4	10	25	39	nat		14	14	14	42
2	8	32	42	nat		15	15	15	45
6	12	24	42	nat		16	16	16	48
1	6	36	43	nat		17	17	17	51
7	14	28	49	nat		18	18	18	54
9	15	25	49	nat		19	19	19	57
4	12	36	52	nat		20	20	20	60
8	16	32	56	nat		21	21	21	63
1	7	49	57	nat		22	22	22	66
12	18	27	57			23	23	23	69
16	20	25	61			24	24	24	72
2	10	50	62	nat		25	25	25	75
3	12	48	63	nat		26	26	26	78
9	18	36	63	nat		27	27	27	81
5	15	45	65	nat		28	28	28	84
4	14	49	67	nat		29	29	29	87
10	20	40	70	nat		30	30	30	90
1	8	64	73	nat		31	31	31	93
18	24	32	74			32	32	32	96
16	24	36	76			33	33	33	99
11	22	44	77	nat
6	18	54	78	nat
8	20	50	78	nat
9	21	49	79	nat
4	16	64	84	nat
12	24	48	84	nat
2	12	72	86	nat
1	9	81	91	nat
7	21	63	91	nat
13	26	52	91	nat
25	30	36	91
3	15	75	93	nat
16	28	49	93	nat
20	30	45	95
9	24	64	97	nat
14	28	56	98	nat
18	30	50	98	nat
```

Edit : nat = not a triangle

Edited: 12 Mar 2012, 11:46 a.m. after one or more responses were posted

Quote:
I'm getting 83 too (see above). They all satisfy the condition as stated.

Yes, you're on the right track. Your list satisfies every condition except one very important one... See my question above.

Quote:
Yes, all the 83 I count satisfy b2=ac condition.
I also get 50 by considering only a < b < c.

Ok, I assume you also get the 1,2,4 shape as above? What is the largest angle in a triangle with sides 1,2 and 4?

Quote:
I'm getting 83 too (see above). They all satisfy the condition as stated.

- Pauli

As has already been said, the original statement explicitly says we're looking for triangles, not mere triplets of integers, and in every triangle you have that none of the sides can equal or exceed the sum of the others.

So, if you've got alleged triangle sides A, B, and C, then they must satisfy:

A+B > C

A+C > B

B+C > A

else you can't construct a triangle having those sides. If in doubt, just try to draw a triangle having sides 1, 2, and 4, and see if you succeed. You won't.

With that restriction in place you get exactly 42 integer triplets which can indeed form a triangle whose perimeter doesn't exceed 100.

Best regards from V.

Edited to correct a typo, as kindly pointed out below

```
```

Edited: 11 Mar 2012, 7:38 p.m. after one or more responses were posted

OK.

I found 83 triplets (a,b,c), but I have not check at all if there are triangles !

In fact (1,2,4) is not a triangle, I found no way to draw it. One of the edge is too short.

That what I miss in the problem; a triplet (a,b,c) is not necessery a triangle.

Thanks.

I have to found a way to check for triangularity and to adjust my algorithm.

'----------------------------------

```00 { 59-Byte Prgm }
01>LBL "FTRI"
02 STO 01              @  Store perimeter in R01
03 3
04 ÷
05 IP                  @ Initate a to INT(P/3)
06 0                   @ Initiate counter R00 to zero
07 STO 00
08 Rv
09>LBL 00              @ ---------- Main loop (over a in stack x:)
10 RCL ST X
11>LBL 01              @ ---------- Loop over b from a to ...
12 RCL ST X
13 RCL+ ST Z           @    place a+b in stack
14 RCL ST Y
15 X^2
16 RCL÷ ST T           @   estimate c = b.b/a
17 X>Y?
18 GTO 02              @   test triangle if ( c > a+b ) exit
19 +
20 RCL 01
21 X<Y?
22 GTO 02              @  test perimeter if ( a+b+c > P ) exit
23 +
24 FP
25 X=0?                @  test c integer if c interger count triangle
26 ISG 00
27 NEG
28 Rv
29 1                   @ increase b
30 +
31 GTO 01
32>LBL 02              @ -------  exit
33 R^
34 DSE ST X
35 GTO 00              @     loop while a>0
36 RCL 00              @ ------ recall counter
37 END
```

Register:

R00 store n (number of geometric triangles)

R01 store perimeter P

usage:

100 EXQ FTRI ---> 42

1000 EXQ FTRI ---> 532

10000 EXQ FTRI ----> 6427

Edited: 11 Mar 2012, 9:01 p.m.

Quote:
So, if you've got alleged triangle sides A, B, and C, then they must satisfy:

A+B < C

A+C < B

B+C < A

If you change all 3 '<' to '>', then you're right. ;-)

Yes, of course, thank you.

Unless I incur in some other typo, I get (in reasonable times):

```   F(100)   = 42
F(1000)  = 532
F(10000) = 6427
```
Best regards from V.

Valentin,
Greetings, What platform are you using? The trusty 71b?

I knew I was missing something obvious, I feel silly now :-)

I've not been thinking well recently :-(

- Pauli

An alternative to the usual triangle inequality arises from Heron's formula, which makes it clear that the semiperimeter of an actual triangle is greater than all three side lengths.

Quote:
10000 EXQ FTRI ----> 6427

A very compact and small size solution. (Thank you for commenting the various loops and code) Very Nice!

Now as to the speed... How long does it take your emulator to calculate F(1e6)? What about F(1e8)?

Is there a clever way to speed up your program by a factor of 1000 or more?

Here are some target times for different machines:

```F(N)        41C     32sii       42s    Free42(D)    Result
------------------------------------------------------------
F(100)      9s        2s        6s        0s            42
F(1000)    72s       15s       45s        0s           532
F(1e4)    651s      145s      426s        0s         6,427
F(1e5)      -         -         -         0s        DD,DDD
F(1e6)      -         -         -       1.5s       CCC,CCC
F(1e7)      -         -         -         5s     B,BBB,BBB
F(1e8)      -         -         -        60s   AAA,AAA,AAA
------------------------------------------------------------
```

Quote:
Valentin,
Greetings, What platform are you using? The trusty 71b?

Hi, Allen.

Yes, I'm using the ole trusty HP-71B in its J-F Garnier's Emu71 incarnation, as usual. Using the naïve approach you can't do better than this:

```    >LIST
10 DESTROY ALL @ INPUT K @ L=K DIV 3 @ N=L @ F=(1+SQR(5))/2
20 SETTIME 0 @ FOR A=1 TO L @ M=(K-A)*A @ FOR B=A+1 TO MIN(A*F,L)
30 IF MOD(B*B,A) THEN 40 ELSE IF B*(A+B)>M THEN 50 ELSE N=N+1
40 NEXT B
50 NEXT A @ DISP N,TIME
>RUN
? 100
42    .02
>RUN
? 1000
532    .52
>RUN
? 10000
6427   51.57
```

which uses two nested loops and minimizes the number of arithmetic operations and comparisons in the inner one to speed execution up.

Nevertheless, using two loops implies increasing the running time by two orders of magnitude (x100) when the input increases by just one (x10), as was bound to be the case and as the timings above clearly demonstrate thus making this naïve approach utterly hopeless for large inputs.

What's needed is a subtler approach that keeps the running time approximately linear on the input, not quadratic as shown here, but I won't spoil the fun for now ... XD

Nice challenge, Allen, thanks for it and best regards from V.

```
```

Edited: 12 Mar 2012, 8:33 a.m.

Quote:
F=(1+SQR(5))/2

I came up with the constraint A > (SQR(5)-1)/2 * B, using a pencil and paper approach. Here are some thoughts:

1. Loop over B (from 1 to N/3, maybe one more and check A+B+C)
2. Add one to the count for the case A=B=C.
3. Loop over A from CEIL((SQR(5)-1)/2 * B) to B-1 and check if B*B/A is an integer.
The last step can be omitted if B is prime or a power of a prime if this can be easily determined.

Edited: 12 Mar 2012, 10:16 a.m.

I made the same mistake too, but I prefer to think of it as a temporary trip into a non-Euclidean geometry realm :-)

Gerson.

I did almost include a comment about the problem not specifying the triangles live in a Euclidean space. Riemannian spaces are so much more fun and correspond better to the word we live in.

- Pauli

Quote:
Here are some thoughts:
1. Loop over B...
2. Loop over A...

Yes, this is in the right direction, Keep in mind, when dealing with
Big O Notation reducing the limits may be slightly faster, but not exponentially faster, since you have to loop over BOTH A and B to solve F(N).

As N grows, your solve time will still grow quadratically O(N^2)- pronounced "Order N squared".

I have studied this problem for over a month- and I am convinced you could teach a large part of an undergraduate computer science course based on solving this problem. In the initial Problem 370, the writers ask for F(2.5e13). The 14-digit answer is impractical on a pocket calculator. But smaller limits like F(1e8) are well within "reasonable" times for emulated calculators, and F(10,000) within reasonable times on a real calculator.

After a month of study, I've seen/found/broken four different approaches to this problem:

1. Naive approach O(N^2) with two loops (easy to write but slow) - like Valentin's nice 71b solution
2. Reduced O(N^2) approach with two loops, with more mathematical limits and assumptions- See Gerson and C.Ret's short and well posed solution- I like them both for different reasons.
3. An O(N) approach (I will post shortly)
4. The most common approach posted in the Euler Discussion forum, which I will not discuss here.

Except for one program, all of the posted approaches in the Project Euler discussion forum for Problem 370 are in category 4 (due to the large size of 2e13). I believe category 4 approach is quite impractical on a pocked calculator for reasons I won't discuss here. But there is at least 1 category 3 approach that is worth discussing in more detail....

Edited: 12 Mar 2012, 9:09 p.m.

Quote:
Apart from the trivial cases of a=b=c this program would seem to produce (I converted it to C to play with so could be wrong here) these nine extras:

``` 25  30  36
20  30  45
18  24  32
16  24  36
16  20  25
12  18  27
9  12  16
8  12  18
4   6   9
```

Your conversion to C is perfect. Here is the BASIC program I wrote to test my second algorithm. The RPN code is just a quick & dirt conversion from it.

```10 CLS
15 INPUT L: PRINT
20 N = INT(L / 3)
25 B = N
30 A = B - 1
35   C = B * B / A
40   IF INT(C) > (A + B) THEN 65
45   IF C - INT(C) <> 0 THEN 55
50     IF C + A + B <= L THEN N = N + 1:  PRINT USING "###"; A; B; C
55     A = A - 1
60   GOTO 35
65   B = B - 1
70 IF B > 2 THEN 30
75 PRINT : PRINT N
? 100
25 30 36
20 30 45
18 24 32
16 24 36
16 20 25
12 18 27
9 12 16
8 12 18
4  6  9
42
```

Would you mind posting your C program? Yesterday I tried to convert it to C by I quit after a while (literally after a while :-). Looks like my BASIC program was not well structured enough. I changed it a bit later but still had trouble with the inner loop when trying to convert it to C:

```10 CLS
15 INPUT L
20 N = INT(L / 3)
25 FOR B = N TO 2 STEP -1
30   A = B - 1
35     C = B * B / A
40     IF INT(C) > (A + B) THEN 65
45     IF C - INT(C) = 0 THEN IF C + A + B <= L THEN N = N + 1
55     A = A - 1
60   GOTO 35
65   '
70 NEXT B
75 PRINT N
```

Thanks,

Gerson.

Seems the code had disappeared, but here is a second stab which is pretty close to the original.

- Pauli

```#include <stdio.h>
const int p = 100;
int main() {
int a, b;
const int n = p / 3;
int count = n;
for (b=n; b>0; b--)
for (a=b-1; a>0; a--) {
const int c = b*b/a;
if (b*b == a*c && c <= a+b && a+b+c <= p)
count++;
}
printf("%d\n", count);
return 0;
}
```

Thank you!

F(100 000) = 75243 in 10 seconds. F(1000 0000) will take a while more * :-)

Gerson.

* this would require int_64, I think.

Edited: 12 Mar 2012, 10:27 p.m.

Yes it will. (10,000,000 / 3)2 exceeds the size of a 32 bit integer. Even 1,000,000 will be too large.

Given that the process is O(n2), I'd expect 1,000,000 to take 100 times as long and 10,000,000 about 10,000 times as long (very roughly).

I've been trying to think of a way of iterating over c and determining possible values of a & b to scan over. No luck thus far. We know that a.x2 = b.x = c for some values of x.

- Pauli

Hi, Allen:

Quote:
After a month of study, I've seen/found/broken four different approaches to this problem:

After a couple' hours of study I came up with a O(n) solution (or O(n)1+e), or O(n*ln(n)), difficult to assess exactly ...) which can find the solution for reasonably large values of N in reasonable times on Emu71, about 1,000 times faster in a compiled high-level language.

I suggest we all post and comment our solutions next Friday so as not to spoil the fun for people working right now on their very own approach. Deal ? :D

Best regards from V.

```
```

Quote:
I came up with the constraint A > (SQR(5)-1)/2 * B, using a pencil and paper approach.

Since

```b^2 = a*c
```
and
```c <= a + b
```
then
```b^2 <= a*(a + b)
b^2 <= a^2 + a*b
b^2/a^2 <= 1 + b/a
(b/a)^2 - b/a - 1 <= 0
```
Solving for b/a:
```b/a <= (sqrt(5) + 1)/2
a/b > (sqrt(5) - 1)/2
a > b*(sqrt(5) - 1)/2
```

I've modified my algorithm to use this constraint. I thought there was no need to check the sum of the sides of the triangles because that's what this constraint was supposed to, but it had to stay. Here are the corresponding C program and RPN code:

```#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main()
{
int a, b, len, n;
scanf("%d", &len);
double c, k;
n = len/3;
k = (sqrt(5) - 1) / 2;
for (b = n; b >= 2; b--)
{
a = b - 1;
while (a > b * k )
{
c = 1.0 * b * b / a;
if (c - int(c) == 0)
if (a + b + c <= len)
n++;
a--;
}
}
printf("n = %d\n",n);
system("pause");
return 0;
}
00 { 85-Byte Prgm }
01>LBL "GT"
02 STO 00                       R00 <- len
03 3
04 BASE÷
05 1E-3
06 +
07 STO 01                       R01 <- b
08 IP
09 STO 02                       R02 <- n = int(len/3)
10 5
11 SQRT
12 1
13 -
14 2
15 ÷
16 STO 04                       R04 <- k = (sqrt(5) - 1)/2
17>LBL 01                       for b = n down to 2
18 RCL 01
19 1
20 BASE-
21 STO 03                       R03 <- a = b - 1
22>LBL 02
23 RCL 01
24 IP
25 RCL× 04
26 RCL 03
27 X<=Y?                         a <= b*k ?  (while a > b*k)
28 GTO 04
29 RCL 01
30 IP
31 X^2
32 RCL÷ 03                      c = b^2/a
33 FP
34 X!=0?                        frac(c) != 0 ?
35 GTO 03
36 LASTX                          goto 03
37 RCL 01                        else
38 BASE+
39 RCL+ 03
40 RCL 00
41 X<=Y?                         c + b + a >= len?
42 GTO 03                          goto 03
43 1                              else
44 STO+ 02                         n = n + 1
45>LBL 03
46 1
47 STO- 03                      a = a - 1
48 GTO 02                       endwhile
49>LBL 04
50 DSE 01
51 GTO 01                       endfor
52 RCL 02                       disp n
53 .END.
100 XEQ GT  -->    42     (1 min 10 sec, real HP-42S)
1000 XEQ GT  -->   532     (~ 0.5 seconds, Free42 Decimal)
10000 XEQ GT  -->  6427     ( 20.8 seconds, Free42 Decimal)
100000 XEQ GT  --> 75243     (~ 41 minutes,  Free42 Decimal)
```
Still an O(n2) process. I've been busy these days and haven't given it much thought, if this can be an excuse. Anyway, no progress after three days. Perhaps I have to upgrade my hardware (brainware) :-)

Gerson.

Edited to fix a typo as per Valentin's observation below.

Edited: 14 Mar 2012, 12:14 p.m. after one or more responses were posted

This is faster, but still an O(n2) process:

```#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int main()
{
int a, b, len, n;
scanf("%d", &len);
double c, k;
n = len/3;
k = (sqrt(5) - 1) / 2;
for (b = n; b >= 2; b--)
{
a = b - 1;
while (a > b * k )
{
c = 1.0 * b * b / a;
if (c - int(c) == 0)
if (a + b + c <= len)
n++;
a--;
}
}
printf("n = %d\n",n);
system("pause");
return 0;
}
100000
n = 75243     (7 seconds @ 1.86 GHz)
```

Gerson.

Hi, Gerson:

Sorry but methinks you've got a typo (unless it's a program bug, hopefully not):

Quote:
``` 10000 XEQ GT  -->  6437     ( 20.8 seconds, Free42 Decimal)
```

The correct value is 6427, not 6437.

Best regards from V.

```
```

My paper and pencil approach was just what you found out: find the zero of a+b>c with ac=b^2. This doesn't make sure that a+b+c is bound in any way, that's just another inequality to check for.

That was indeed a typo. The program does return 6427. I will fix it presently. Thank you for pointing it out.

Best regards,

Gerson.

Edited: 14 Mar 2012, 12:14 p.m.

That's right !

I am really glad to learn that my first approach doesn’t below to the first naïve category and is consider as ‘reduce naïve’.
Because, in my point of view, my first code is really a “brute force scan”, since very few optimization in the search strategy is done.
Two loops imbedded.
By luck the outer one scan step by step on a which is fortunately the shortest side of the triangle, make it walk trough it as short as possible because limited to (P/3) ! How Lucky I am here!
On the other hand, the inner loop (which will weight much more of the time to find the solution) is far for a good optimization. Despite the two exit tests which appropriately short-cut the long b walk through, the worst in my code is that at b is scan step by step (only increase by one). Here I am bad.
This bad strategy is perhaps not a shame for small perimeter, but will be a disaster for larger perimeter.

As a lot of Euler problems, here a more tricky strategy have to be take into account, to avoid this very long systematic scans of “brute force like” algorithm.

There must be a way to scan over a only, not over b !

What are the constraints:

(1) Sorted Sides: a <= b <= c

(2) Being a triangle c <= a + b

Since other inequalities are both verified:

Due to a < b we get a <= b + n whatever is the natural integer n.
Due to initial condition (1) b <= c so we get b <= c + n for any natural integer n, thus b <= c +a.
As long as initial condition (1) is respected, we only have to care about c <= a + b

(3) Limited perimeter a+b+c <= P
As Marcus Von Cube pointed it out, this constraint has to be tested. But generally it is less constringent than (2).
(4) Specific condition: b2 = a.c
Here is the key to reduce

Suppose that b decompose into the following product of primary factors : b1, b2, b3

With b = b1. b2. b3 then (4) implies that a and c are both a subset of the primary factors of b.

To illustrate that, investigate all the triplets (a,b,c) that verify (4). Only a fraction of this list of triplets corresponds to constraints (1) (2) and (3):

```Looking for triangle P<=100
b = 2 x 2 x 3 = 12 scanning all (a,c) pairs built up table 1:
TABLE 1:
========
a    b    c   P=a+b+c Violate constraints
----- ---- ---- --------- ------------------------------------------
1   12  144     157   (2)(3)      Not a triangle + P Too large
2   12   72      86   (2)         Not a triangle
3   12   48      63   (2)         Not a triangle
4   12   36      52   (2)         Not a triangle
6   12   24      42   (2)         Not a triangle
8   12   18      38               Ok n°1
9   12   16      37               Ok n°2
12   12   12      36               Ok n°3
16   12    9      37   (1)         a b c not sorted
18   12    8      38   (1)         a b c not sorted
24   12    6      42   (1)         a b c not sorted
36   12    4      52   (1)         a b c not sorted
48   12    3      63   (1)         a b c not sorted
72   12    2      86   (1)         a b c not sorted
144   12    1     157   (1)(3)      a b c not sorted + Too large P
----- ---- ---- --------- ------------------------------------------
```
We can observe :
- the symmetry role of a and c, scanning over a or c is equivalent.
- half of the triplet have a>b (i.e. c<b), scan can be reduce considerably (for example by only scanning over a from 1 to b)
- not every avalue have to be scan.
In fact, only value of a composed using a subset of the primary factor from b2 :
```Table 2 :
=========
b  =  2 x 2 x 3            b.b =  2 x 2 x 2 x 2 x 3 x 3
------------------------- ------------------------------
a  =  1                      c =  2 x 2 x 2 x 2 x 3 x 3 = 144
a  =  2                      c =      2 x 2 x 2 x 3 x 3 =  72
a  =  3                      c =  2 x 2 x 2 x 2 x 3     =  48
a  =  2 x 2                  c =          2 x 2 x 3 x 3 =  36
a  =  2 x 3                  c =      2 x 2 x 2 x 3     =  24
a  =  2 x 2 x 2              c =              2 x 3 x 3 =  18
a  =  3 x 3                  c =  2 x 2 x 2 x 2         =  16
a  =  2 x 2 x 3              c =          2 x 2 x 3     =  12
...
```
This prove that scanning step by step over a or b is not the trickiest way to look for such triangle.

This suggest that an algorithm based on combination of a prime factors decomposition of bcan be more adapt to this Week-End challenge.

On the other hands, having to scan over the smallest side a will be much simpler and allow the uses of short-cut tests. I am on the way thinking that investigating how to built the correct b side from a given a under constrains (1) to (4) will be an elegant and surely efficient method.
Of course, the over-head need to analyze a and built up b will only be benefic for scanning over large Perimeter and large side. But the great interest would be that no more inner loop over b will be needed.

To investigate how to construct b from a given a, please consider a first example:
Looking for geometric triangle responding to constraints (1) to (4) with perimeter less or egal to P = 100.
Considering a prime :

```Table 3:   prime  a = 7  , P <=100
=======
a    b    c   P=a+b+c Violate constraints
----- ---- ---- --------- ------------------------------------------
7    7    7      21   none                OK
7   14   28      49   (2)                 not a triangle (c>a+b)
7   21   63      91   (2)                 not a triangle
7   28  112     147   (2)(3)              not a triangle & too large
7   35  175     217   (2)(3)              not a triangle & too large
and so on...
----- ---- ---- --------- ------------------------------------------
```

Only first triplet (7,7,7) correspond to the researched geometric triangle.
As soon as the second line, triplets violate one or more constrains, so I may have stop the scan as soon.
Putting more values in the table shows how b value are build.
As a is prim, there is no combination of factors is possible, and we observe that b is multiple of a
The triplets are of the general formulae (a , k.a , a.k2) with k=1, 2, 3, ...

Now, what if a is a simple composite: a = 8

```Table 4:   simple composite a = 8 , P <=100
=======
a    b    c   P=a+b+c Violate constraints
----- ---- ---- --------- ------------------------------------------
8    8    8      24   none                OK
8   12   18      38   none                OK
8   16   32      56   (2)                 not a triangle
8   20   50      78   (2)                 not a triangle
8   24   72     104   (2)(3)              not a triangle & too large
8   28   98     134   (2)(3)              not a triangle & too large
and so on...
----- ---- ---- --------- ------------------------------------------
```

Here, two triangles of small side 8 exist (8,8,8) and (8,12,18). Both candidates respect all the constraints.

If we consider the b set , we now observe an arithmetic suite a, a+4, a+2x4, a+3x4,…
At each rwo, we have to increase b by 4. Where come this value from?
Is there any chance that 4 is built from half of the prime factors of a ?
a = 1x2x2x2 = 8 , the delta D to use to built the set of b is D = 2x2 = 4

The triplets are of the form (a, a+D.k, c)., where k = 0, 1, 2, ....

Now, what if a is a[italic] square: [italic]a = 9

```Table 5:   simple composite a = 25 , P <= 100
=======
a    b    c   P=a+b+c Violate constraints
----- ---- ---- --------- ------------------------------------------
25   25   25      75   none                OK
25   30   36      91   none                OK
25   35   49     109   (3)                 perimeter too long
25   40   64     129   (3)                 perimeter too long
25   45   81     151   (2)(3)              not a triangle & too long
25   50  100     175   (2)(3)              not a triangle & too long
25   55  121     201   (2)(3)              not a triangle & too long
and so on...
----- ---- ---- --------- ------------------------------------------
```

At each row, b have to be increase by D = 5 (not 25). So we have to construct D as small as possible with “half” of the factors of a.
Triplets are of the form ( a2, a2+k.a, (a+k) 2).

In fact, D has to be built from each factor of a, but using only half of the multiple factors (i.e; with half to power of each factor).

```Table 6:   simple composite a = 42 , P <= 500
=======
a    b    c   P=a+b+c Violate constraints
----- ---- ---- --------- ------------------------------------------
42   42   25     126   none                OK
42   84  168     294   (2)                 not a triangle
42  126  378     546   (2)(3)              nat & perimeter too long
42  168  672     882   (2)(3)              nat & perimeter too long
----- ---- ---- --------- ------------------------------------------
```

Here a = 1x2x3x7 = 42 , and b = a + 42.k with k = 0, 1, 2, ....
This confirming that D has to be built from each factor of a, but using only half of the multiple factors (i.e; with half to power of each factor).

```Table 7:   Samples of D from factor decomposition of a
=======
a                    D              |   a                    D               | a                       D              |
-----   ------------ ----- -------------|-----   ------------ ----- -------------|-----   ------------ ----- -------------|
1   1x1              1              |   31   1x31            31              |   61   1x61            61              |
2   1x2              2              |   32   1x2x2x2x2x2      8  = 2x2x2     |   62   1x2x31          62  = 2x31      |
3   1x3              3              |   33   1x33            33              |   63   1x63            63              |
4   1x2x2            2  =  2        |   34   1x2x17          34  = 2x17      |   64   1x2x2x2x2x2x2    8  = 2x2x2     |
5   1x5              5              |   35   1x5x7           35  = 5x7       |
6   1x2x3            6  = 2x3       |   36   1x2x13          36  = 2x2x2     |
7   1x7              7              |   37   1x37            37              |
8   1x2x2x2          4  = 2x2       |   38   1x2x19          38  = 2x19      |
9   1x3x3            3  =  3        |   39   1x3x13          39  = 3x13      |
10   1x2x5           10  = 2x5       |   40   1x2x2x2x5       20  = 2x2x5     |
11   1x11            11              |   41   1x41            41              |
12   1x2x2x3          6  = 2x3       |   42   1x2x3x7         42  = 2x3x7     |
13   1x13            13              |   43   1x42            43              |
14   1x2x7           14  = 2x7       |   44   1x2x2x11        22  = 2x11      |
15   1x3x5           15  = 3x5       |   45   1x3x3x5         15  = 3x5       |
16   1x2x2x2x2        4  = 2x2       |   46   1x2x23          46  = 2x23      |
17   1x17            17              |   47   1x47            47              |
18   1x2x3x3          6  = 2x3       |   48   1x2x2x2x2x3     12  = 2x2x3     |
19   1x19            19              |   49   1x7x7            7  =  7        |
20   1x2x2x5         10  = 2x5       |   50   1x2x5x5         10  = 2x25      |
21   1x3x7           21  = 3x7       |   51   1x3x17          51  = 3x17      |
22   1x2x11          22  = 2x11      |   52   1x2x2x13        26  = 2x13      |
23   1x23            23              |   53   1x53            53              |
24   1x2x2x2x3       12  = 2x2x3     |   54   1x2x3x3x3       18  = 2x3x3     |
25   1x5x5            5  =  5        |   55   1x5x11          55  = 5x11      |
26   1x2x13          26  = 2x13      |   56   1x2x2x2x7       28  = 2x2x7     |
27   1x3x3x3          9  = 3x3       |   57   1x3x19          57  = 3x19      |
28   1x2x2x7         14  = 2x7       |   58   1x2x29          58  = 2x29      |
29   1x29            29              |   59   1x59            59              |
30   1x2x3x5         30  = 2x3x5     |   60   1x2x2x3x5       30  = 2x3x5     |
-----   ------------ ----- -------------|-----   ------------ ----- -------------|
```

The larger is D, the more rapid is the scan. As we can see, the process (involving the two embedded loops) has a great chance to be speed-up, since most of the a value leads to D = a. Especially when a is prime.

The inner loop (the one over b) can be suppress if we found a way of forecasting at which k we have to stop because the triplet violate one (or more) constraints.

To test for constraint (2) : ( a, b, c) being a triangle, we have to test c >= a + b

From the given a value, we may get D (see table 7 above) or find a way to built D from a factors.
With such a D, we know that the triplet (a, a+k.D, c) respect (4) b2 = a.c when
b = a + k.D , therefore c = b2 / a

To test for constrinat (2), we have to express c >= a + b as a function of variable k and parameters a and D.

The inequality c >= a + b can be transform into k2.D2 + k.D.(2-a) - a2 <= 0. Resolving the corresponding quadratic equation, allow to indicate the domain of validity:

The variable have to be set from 0 up to k2 = -(a/2D).(1-V5)

To test for constraint (3) : a+b+c <=P , we may found the limit for k the same way, by expression the inequality as a function of variable k and parameters a, P and D
a+b+c <=P is transform into k2.D2 + k.D.(2+a) + 3a2-aP <= 0.
Solving the corresponding quadratic equation for k indicate that k may varie from 0 upto k3 = -(a/2D).(3-V(4P/a - 3))

This give you the way to code for a more efficient algorithm using only one loop over a :

```Input   P  upper limit of perimeter
Initiate Sum = 0
For each integer a from 1 to P/3
Determine D  ( 1 <= D <= a ) from the prime factor decomposition of a.
In most of the cases D will be equal to a
or less than a in specific cases
IF P > (3+SQR(5))*a
THEN   ' Limited by triangle (2)
kl = -(a/2D)*(1-SQR(5))
ELSE   ' Limited by perimeter (3)
kl = -(a/2D)*(3-SQR(4P/a - 3))
End if
There is n = INT( 1 + kl ) geometric triangles of small side [italic]a[/italic] since k = 0,1, 2, ..., kl
Sum = Sum + n         Add integer n into the Sum
At end, display Sum that contains number of researched triangle of perimeter less or equal to P.
```

In conclusion, I leave you a few to compose you own code on your favorite calculator or system.
I curious to know how much time this new algorithm loose on shot perimeter (P = 100, 1000), and what it is win for larger size (P> 10^4).

Edited: 15 Mar 2012, 7:26 p.m. after one or more responses were posted

Quote:
I am really glad to learn that my first approach doesn’t below to the first naïve category and is consider as ‘reduce naïve’. Because, in my point of view, my first code is really a “brute force scan”, since very few optimization in the search strategy is done.

That is exactly my feelings, regarding mine. I haven't timed your first approach, but mine isn't any better than Valentin's, quite the contrary:

```10 DESTROY ALL @ INPUT L @ N= L DIV 3 @ K=(SQR(5)-1)/2
20 SETTIME 0 @ FOR B = N TO 2 STEP -1 @ A = B - 1
30 IF A < B * K THEN 90
40 C = B * B / A
50 IF FP(C) <>  0 THEN 70
60 IF A + B + C <= L THEN N = N + 1
70 A = A - 1
80 GOTO 30
90 NEXT B @ PRINT N,TIME
>RUN
? 100
42                   .05
>RUN
? 1000
532                  2.51
>RUN
? 10000
6427                 255.22
```
Valentin's times in my notebook running EMU71 @ 1.86 GHz are:

```>RUN
? 100
42                   .05
>RUN
? 1000
532                  1.2
>RUN
? 10000
6427                 114.67
```

Gerson.

You may try this last version :

```(Microsoft Q-BASIC)
DEFLNG A-Z         '  ---- All variables long integer except k and t (single precision)
DIM k AS SINGLE
DIM t AS SINGLE
PRINT
INPUT "ENTER perimeter"; P
IF P < 1 THEN END
'
'   Initiate Sum and timer
'
Sum = 0
t0 = TIMER
'
'  --------------------- MAIN LOOP OVER a
'
FOR a = 1 TO P / 3
'
' ------------------------ Compute step D
'
x = a: i = 2: j = 1: D = 1
'
'  - - - - - - - - - - Isprime? Test Loop
DO UNTIL i * i > x
n = 0
DO WHILE x MOD i = 0: n = n + 1: x = x / i: LOOP
IF n > 0 THEN D = D * i ^ INT(.5 + n / 2)
i = i + j: j = 2
LOOP
D = D * x
'
' ----------------------- test Limit for k
'
IF P < (3 + SQR(5)) * a THEN
k = SQR(4 * P / a - 3) - 3   ' Limited by perimeter
ELSE
k = SQR(5) - 1               ' Limited by triangle
END IF
'
' ----------------------- Add current number to sum
'
Sum = Sum + 1 + INT(a / 2 / D * k)
NEXT a
'
'   ------------------------- Display result and time
'
PRINT : PRINT
PRINT "There is "; Sum; "triangles of perimeter <="; P; " with b²=a.c"
PRINT "("; INT((TIMER - t0) * 10) / 10; "s )"
PRINT
```

Approximetely translate into HP-71 BASIC : not sure about correct syntax of MOD CEIL and power (^).

``` 10 DESTROY ALL @ INPUT P @ SETTIME 0
REM -------------  MAIN LOOP OVER A
20 S = 0 @ FOR A = 1 TO P DIV 3
REM -------------  DETERMINE STEP D
30    X = A @ I = 2 @ D = 1 @ J = 1
REM - - - - - - - ISPRIME TEST LOOP
40    N = 0 @ IF I * I > X THEN 90
50       IF X MOD I = 0 THEN N = N + 1 @ X = X DIV I @ GOTO 50
60       IF N > 0 THEN D = D * I ^ CEIL( N / 2 )
70       I = I + J @ J = 2
80    GOTO 40
REM -------------  COUNT SOLUTIONS
90    D = D * X @ K = 1 - SQRT 5
100    IF P < (3 + SQRT 5) * A THEN K = 3 - SQRT(4 * P / A - 3)
110    S = S + 1 + INT( - A / 2 / D * K)
120 NEXT A @ PRINT S,TIME
```

Edit : Replace N by S at line 120 in HP71 listing.

Edited: 15 Mar 2012, 6:53 p.m. after one or more responses were posted

It looks like there is an error in your conversion from QBASIC to HP-71 BASIC or I have introduced one myself:

```10 DESTROY ALL @ INPUT P @ SETTIME 0
20 S=0 @ FOR A=1 TO P DIV 3
23 REM
25 REM -------------  determine step d
27 REM
30 X=A @ I=2 @ D=1 @ J=1
33 REM
35 REM - - - - - - - isprime test loop
37 REM
40 N=0 @ IF I*I>X THEN 90
45 N=0
50 IF MOD(X,I)=0 THEN N=N+1 @ X=X DIV I @ GOTO 50
60 IF N>0 THEN D=D*I^CEIL(N/2)
70 I=I+J @ J=2
80 GOTO 40
83 REM
85 REM -------------  count solutions
87 REM
90 D=D*X @ K=1-SQR(5)
100 IF P<(3+SQR(5))*A THEN K=3-SQR(4*P/A-3)
110 S=S+1+INT(-A/2/D*K)
120 NEXT A @ PRINT N,TIME
>RUN
? 1000
0                    .32
>RUN
? 10000
0                    5.1
>RUN
? 100000
0                    102.21
```

The QBASIC times are:

```    F(100000) -->   1.4 s
F(1000000) -->  30.6 s
F(10000000) --> 567.8 s
```

Yes, I made a mistake in the last line where N stand for S.
Sorry.

```10 DESTROY ALL @ INPUT P @ SETTIME 0 @ S=0 @ FOR A=1 TO P DIV 3 @  X=A @ I=2 @ D=1 @ J=1
40 N=0 @ IF I*I>X THEN 90
50 IF MOD(X,I)=0 THEN N=N+1 @ X=X DIV I @ GOTO 50
60 IF N>0 THEN D=D*I^CEIL(N/2)
70 I=I+J @ J=2 @ GOTO 40
90 D=D*X @ K=SQR(5)-1 @ IF P<(3+SQR(5))*A THEN K=SQR(4*P/A-3)-3
110 S=S+1+INT(A/2/D*K) @ NEXT A @ PRINT S,TIME
```

Here the same code for "old" RPL aka HP-28C/S.
"New PRL" may take advantage of ISPRIME? and NEXT-PRIME instructions !

```« -> P
« 0                                        @ Initiate Sum=0
1 P 3 / FOR a
1 SF
2                                    @ Initiate D = 2
a                                    @          x = a
2                                    @          i = 2
WHILE DUP2 SQ >=                     @ while i²<=x
REPEAT
0 3 ROLLD                         @    Initate  n = 0
WHILE DUP2 MOD NOT                @    while x MOD i = 0
REPEAT
SWAP OVER / SWAP               @       x = x/i
ROT 1 + 3 ROLLD                @       n = n+1
END
4 ROLL OVER 5 ROLL 2 /
CEIL ^ * 3 ROLLD                  @       D = D*i^CEIL(n/2)
1 +                               @   Inc i
IF 1 FC?C THEN 1 + END            @   Inc i when i>3
END
DROP *                               @   D = D*x
a SWAP /
IF P a 3 5 SQRT + * <
THEN 4 P * a / 3 - SQRT 3 -       @ limited by perimeter
ELSE 5 SQRT 1 -                   @ limited by triangle side
END
* IP 1 + +                           @   Sum = Sum + 1 + INT(...)
NEXT
»
»  'B2AC' STO
```

Execution Times HP28S
100 B2AC ---> 42 (12s.)
1000 B2AC ---> 532 (3min45s)

Edited: 15 Mar 2012, 8:04 p.m. after one or more responses were posted

Here are your times up to 10^5:

```  Emu71: HP-71B & HP-IL system emulator         J-F GARNIER 1996, 2006
>list
10 DESTROY ALL @ FOR E=2 TO 5 @ P=10^E @ SETTIME 0 @ S=0
20 FOR A=1 TO P DIV 3 @ X=A @ I=2 @ D=1 @ J=1
40 N=0 @ IF I*I>X THEN 90
50 IF MOD(X,I)=0 THEN N=N+1 @ X=X DIV I @ GOTO 50
60 IF N>0 THEN D=D*I^CEIL(N/2)
70 I=I+J @ J=2 @ GOTO 40
90 D=D*X @ K=SQR(5)-1 @ IF P<(3+SQR(5))*A THEN K=SQR(4*P/A-3)-3
110 S=S+1+INT(A/2/D*K) @ NEXT A @ PRINT P,S,TIME @ NEXT E
>run
100                  42                   .05
1000                 532                  .3
10000                6427                 4.81
100000               75243                94.45
```

Thank you.

It's look like this version is not as fast as Valentin Albillo's code (see post #023) who is using machine ROM PRIM function.

My code loose much of his time to test for primility, where the PRIM function spare a lot of time.

As as explain above, most od the time, D equal a. In most of the case my code test all factors of a where a simple D = a will do the job.

Following a transcription of the above code into HP-32S RPN program :

```00 { 141-Byte Prgm } :                     :
01>LBL "B2AC"        : 31 1                : 61 SQRT
02 STO 01            : 32 +                : 62 +
03 3                 : 33 GTO 02           : 63 ×
04 ÷                 : 34>LBL 03<          : 64 X<=Y?
05 IP                : 35 Rv               : 65 GTO 06
06 STO 02            : 36 X=0?             : 66 Rv
07 0                 : 37 GTO 04           : 67 4
08 STO 00            : 38 2                : 68 ×
09>LBL 00<-----------: 39 1/X              : 69 RCL÷ 02
10 2                 : 40 STO× ST Y        : 70 3
11 STO 03            : 41 +                : 71 STO- ST Y
12 RCL 02            : 42 IP               : 72 X<>Y
13 2                 : 43 RCL ST Y         : 73 GTO 07
14 SF 00             : 44 X<>Y             : 74>LBL 06<-------
15>LBL 01            : 45 Y^X              : 75 1
16 X^2               : 46 STO× 03          : 76 5
17 X>Y?              : 47>LBL 04<          : 77>LBL 07<-------
18 GTO 05            : 48 CLX              : 78 SQRT
19 X<> ST L          : 49 1                : 79 X<>Y
20 0                 : 50 FC?C 00          : 80 -
21>LBL 02<           : 51 STO+ ST Y        : 81 RCL× 02
22 RCL ST Z          : 52 +                : 82 RCL÷ 03
23 RCL÷ ST Z         : 53 GTO 01           : 83 IP
24 FP                : 54>LBL 05<----------: 84 1
25 X>0?              : 55 Rv               : 85 +
26 GTO 03            : 56 STO× 03          : 86 STO+ 00
27 Rv                : 57 RCL 01           : 87 DSE 02
28 X<>Y              : 58 RCL 02           : 88 GTO 00
29 STO÷ ST Z         : 59 3                : 89 RCL 01
30 X<>Y              : 60 5                : 90 RCL 00
:                     : 91 END
REgisters
R00 : Sum ,number of triangles
R01 :  P  ,perimeter
R02 :  a  ,small side
R03 :  D  ,
>LBL 00 : Main loop (For a = INT(P/3) downto 1
Factors analysis of a
>LBL 01 : Look for factors of a
>LBL 02 : search for multiple factor of a
>LBL 03 : Built D by adding each factor of a (but half of each multiple factor number).
>LBL 04 : next factor f=2,3,5,...
Compute number of triangle of small side a
>LBL 05 : Nbr of triangle limited by perimeter a+b+c < P
>LBL 06 : Nbr of triangle limited by side c < a+b
>LBL 07 : Add number of traingle of side a to global sum (R00)
```

Edited: 16 Mar 2012, 5:56 p.m.

Hi, all

We're already on Friday so time to post some non-naïve solutions to this challenge.

When I created my above naïve solution I fully realized that having two nested loops, though simple and immediately coded, was not the way to go as you'll inevitably get times of order O(n2). Some clever tricks with the nested loop limits did significantly improve the timing but it still was O(n2).

Regrettably I couldn't dedicate more than a few hours to this but I decided to make the most of the available time and last Tuesday I concocted this much improved 6-line, 234-byte solution for the HP-71B which I'll presently comment:

```     1 DESTROY ALL @ INPUT K @ P=4*K @ N=0 @ F=(1+SQR(5))/2 @ SETTIME 0
2 FOR A=1 TO K/3 @ N=N+FNT(FNF(A)) @ NEXT A @ DISP N;TIME
3 DEF FNT(N)=MIN(A*F,A*((SQR(P/A-3)-1)/2)) DIV N-A DIV N+NOT MOD(A,N)
4 DEF FNF(N) @ M=N @ E=0
5 D=PRIM(N) @ IF D<2 THEN FNF=M/(E*(N=E)+(N#E)) @ END ELSE N=N/D
6 IF E=D THEN M=M/E @ E=0 @ GOTO 5 ELSE E=D @ GOTO 5
```

First, some timings:

```         Perimeter      Triangles        Time (sec).   Time increase
------------------------------------------------------------
100              42              0.03           -
1000             532             0.16          5.2x
10000            6427            1.8          11.3x
100000           75243          19.52         10.8x
1000000          861805        219.62         11.3x
10000000         9712598      2589.75         11.8x
```

• As you can see, the main program is just 2-line long and all it does is to initialize some needed constants so that the ensuing loop will go as fast as possible by removing invariant operations within the loop, then the single loop on the shortest side (A) is executed, which simply calls a user-defined function FNT which returns the number of valid triangles for that value of A and adds this up to the tally, eventually returning the total and the timing when the loop is done for. Simple as can be.

• Line 3 is a nifty user-defined single-line function, FNT, that duly returns the number of valid triangles within limits given by the maximum and minimum possible values for side B, which are cleverly computed so that no tests about triangularity and/or perimeter are necessary within the loop at all, this way:

- the A*F term, where F is the golden ratio, ensures that A, B, and C form inded a triangle, as C is guaranteed to be less than A+B.

- the A*((SQR(P/A-3)-1)/2)) term ensures that the sum of A+B+C is guaranteed to be less than or equal to the given perimeter.

• Lines 4-5 constitute another user-defined function, FNF, that given the value of A computes the periodicity with which B*B is divisible by A. This way there's no need for another nested loop that examines each B in turn to check for divisibility but we can simply tell at once just how many values of B do meet the divisibility condition by simply using this computed periodicity as input for the triangle-counting function FNT.

```
```
Thus, this approach does allow us to get over the dreaded O(n2) to achive something much more like O(n*log(n)), as the results above clearly show ( ~ 11.3x increase in time for every 10x increase in the input value).

Note:

FNF uses the PRIM function from the HP-71B JPC ROM for faster execution. If this function were unavailable the functionality can easily be achieved using plain-vanilla HP-71B BASIC code, which would affect speed by a constant factor (~3x slower) but would not change the order O(n*log(n)) so timing would ultimately be much faster than any O(n2) approach for sufficiently large input values of the perimeter.

That's all on my side, thanks to Allen for an interesting challenge.

Best regards from V.

```
```

Valentin, C.Ret, Gerson- thank you for your interesting solutions. I hope I was clear in saying there is nothing wrong with the naive approach- only an observation that it is the easiest solution to develop and program. Your solutions are certainly well thought out and not "naive" as in a bad.

I admit my approach is not the fastest I've seen (or even close) but I propose a method that neither counts a nor b, but only focuses on counting similar triangles with ratio of b/a. In so doing, There are no perimeter checks, no triangle inequality, and no loops over a or b- in fact there is only 1 loop to generate a series of fractions, and 1 check to make sure that fraction b/a does not pass the golden ratio. This allows execution to proceed in O(N) time - requiring roughly N/16 cycles for F(N).

How?

F(100) can be calculated in 7 iterations of fraction loop that creates coprime numerator and denominator pairs called a Farey Sequence (A close cousin of the Stern-Brocot tree). Using Integer Division, these coprime pairs can be used to quickly sum how many similar triangles are less than the perimeter. A more lengthy description can be found in this PDF file: Solving Project Euler.net Problem 370 on an HP calculator. This PDF also includes Bar codes for HP41, .RAW file if you want to use it in an Free42S emulator, commented HP42S Code, and a C program if you wish to compile on a PC.

Some example Run Times:

```F(N)	41C	32sii	42s	Free42(D)  Result	Iterations	Time increase
F(100)	9s	2s	6s	0s	   42	          7	             -
F(1000)	72s	15s	45s	0s	   532	          64	            9.142857143
F(1e4)	651s	145s	426s	0s	   6427	          619	            9.671875
F(1e5)	-	-	-	0s	   75243	  6262	           10.11631664
F(1e6)	-	-	-	1.5s	   861805	  62710	           10.0143724
F(1e7)	-	-	-	5s	   9712598	  626180	    9.985329294
F(1e8)	-	-	-	60s	   108073540	  6262634	   10.00133189
Note: in MATLAB or C the sum of times for F(2) to F(8)) is around 0.35 seconds
```

The time increase approaches 10X for a 10X increase in F(N) so it's almost exactly O(N) all the way up to higher orders.

Here is the commented HP2S code. It uses 10 registers: (1 for the Total, 1 for the order of the sequence, 1 for Phi, the golden ratio, 1 temporary register to simplify calculations, and 6 registers to iterate the Farey Sequence.)

1. Lines 1-22 set up 10 registers for Main loop
2. Lines 24-34 find the smallest similar Triangle with sides a and b
3. Lines 36-62 Creates the next Farey sequence.
4. Lines 63-66 Check to exit if m2/n2> golden ratio

```00 { 72-Byte Prgm }	Comment
01>LBL 01	Lines 1-22 set up 10 registers for Main loop
02 STO 00	p=Perimeter
03 3
04 ÷
05 SQRT	        Set Order for Recursive Farey Series
06 IP
07 STO 01	=floor(sqrt(P/3))
08 SIGN  	save a byte to get 1 on the stack
09 STO 03	mt=1 - Initial Farey upper bound mt/nt=1/1
10 STO 04	nt=1
11 STO 05	m2=current side a=1
12 STO 06	n2=current side b=1
13 STO 07	n1=1
14 5
15 SQRT
16 +
17 2
18 ÷
19 STO 02	Phi Golden Ratio (the upper limit for the ratio of mt/nt)
20 CLX
21 STO 08	m1=0  Initial Farey lower bound m1/n1=0/1
22 STO 09	Total triangles= 0
23>LBL 02	Main loop
24 RCL 00	Lines 24-34 find the smallest similar Triangle with sides a and b
25 RCL 05
26 RCL 06
27 +
28 X^2	        The smallest perimeter is ps=(a+b)^2-ab
29 RCL 05
30 RCL 06
31 ×
32 -
33 ÷
34 IP
35 STO+ 09	Total=total+floor(p/ps) [if ps>p this adds 0- no branching test needed]
36 RCL 07	Lines 36-62 Create the next Farey sequence.
37 RCL 01	  This guarantees the numbers mt and nt are coprime and thus similar triangles are not
38 +	          counted twice.   Example: For the trivial case where the perimeter p=100,
39 RCL 06	  the order 5 Farey series is: 0/1, 1/1, 6/5, 5/4, 4/3, 7/5, 3/2, 8/5…
40 ÷	         Contributing 0+33+1+1+2+0+5+0= 42 respectively.
41 IP	         h=floor((n1+d)/n2)
42 ENTER
43 ENTER
44 RCL 05
45 ×
46 RCL 08
47 -
48 STO 03	mt=h*m2-m1
49 X<>Y
50 RCL 06
51 ×
52 RCL 07
53 -
54 STO 04	nt=h*n2-n1
55 RCL 06
56 STO 07	n1=n2
57 RCL 05
58 STO 08	m1=m2
59 RCL 03
60 STO 05	m2=mt
61 RCL 04
62 STO 06	n2=nt
63 ÷
64 RCL 02
65 X>Y?	        Lines 63-66 Check to exit if m2/n2> golden ratio
66 GTO 02	Next Farey sequence
67 RCL 09	Otherwise recall total (END)
68 .END.
```

For F(2.5e13) the C program will get the solution to Euler Problem 370 eventually, but since it's O(N) this takes a long time. Nonetheless, because of the relatively small size program and minimal data requirements, finding Similar triangles using the Farey Sequence is reasonable for smaller orders of N. For these reasons, I believe this approach is is moderately suited for many values of N using a real or emulated HP calculator.

Many thanks to those who contributed solutions, you have some clever and compact solutions, all of which I studied and learned from.

Edited: 16 Mar 2012, 9:37 p.m.

Congrats, Allen, this is a really beautiful solution!

For fun I coded it up in RPL+ and ended up with this:

```\<< =n
0 =:m1 =sum
1 =:m2 =:n1 =n2
floor(sqrt(n/3)) =d
goldenRatio EVAL =gr
WHILE m2/n2<=gr
REPEAT
floor(n/((m2+n2)*(m2+n2)-m2*n2)) =+sum
floor((n1+d)/n2) =h
h*m2-m1 =mt
h*n2-n1 =nt
n2 =n1 nt =n2
m2 =m1 mt =m2
END
sum
\>>
```

goldenRatio is a built-in continued fraction. It's EVAL'ed in order to yield a real.

Run-times:

ND1 (iPhone): F(1000): <0.1s; F(1e6): 35s

CalcPad (Mac): F(1000): 0.007s; F(1e6): 1.8s

I used your C code to arrive at this in J‌avaScript:

```function (n) { /*as is*/
var sum=0, maxp = n;
var mt=1, nt=1, m2=1, n1=1, n2=1, h, m1=0;
var d = Math.floor(Math.sqrt(maxp/3));
var gr = 1.6180339887498948482;
while (m2/n2 <= gr ) {
sum += Math.floor(maxp/((m2+n2)*(m2+n2)-m2*n2));
h = Math.floor((n1+d)/n2);
mt = h * m2 - m1;
nt = h * n2 - n1;
n1 = n2;
n2 = nt;
m1 = m2;
m2 = mt;
}
return sum;
}
```

Run-times:

ND1 (iPhone): F(1000): 0.003s; F(1e8): 32s

CalcPad (Mac): F(1000): instant; F(1e8): 3.6s

My goal for RPL+ is to eventually achieve roughly half of JavaScript speed.

Edited: 17 Mar 2012, 8:03 a.m.

Nice code as always, C.Ret!

Quote:
Execution Times HP28S 100 B2AC ---> 42 (12s.) 1000 B2AC ---> 532 (3min45s)

I ran it (unchanged) in ND1, yielding the following run-times:
B2AC(100): 0.2s; B2AC(1000): 2.8s

It's not obvious to me how to change the code to use NEXTPRIME and/or ISPRIME?. Can you show me?

The 50g or ND1 also have FACTORS, which can be used to derive your "D". It's straightforward but not terribly elegant, as you'd still need a loop to do an integer divide on every second element in the factors array (which interleaves factors and their exponents).

Hi again, Allen:

Quote:
Valentin, C.Ret, Gerson- thank you for your interesting solutions. I hope I was clear in saying there is nothing wrong with the naive approach- only an observation that it is the easiest solution to develop and program. Your solutions are certainly well thought out and not "naive" as in a bad.

Well, I wouldn't say my 6-line HP-71B solution above is that naïve in any sense good or bad but I get your point. As kinda proof-of-concept, this is the exact same solution directly converted to just 5 lines of interpreted UBASIC code almost verbatim:

```
10  word -3 : point -2 : for I=2 to 8 : K=10^I : P=4*K : S=0 : F=(1+sqrt(5))/2 : clr time
20  for A=1 to K\3 : N=A : M=N : E=0
30  D=prmdiv(N) : if D>1 then N\=D : if E=D then M\=D : E=0 : goto 30 else E=D : goto 30
40  S+=int(min(A*F,A*((sqrt(P/A-3)-1)/2))/M)-A\M+(res=0) : next
50  print I,S,time : next : print "Yasta"
```

Upon running it, the results and execution times are, in the same hardware:

```run
2   42          0:00:00
3   532         0:00:00
4   6427        0:00:00
5   75243       0:00:00
6   861805      0:00:01
7   9712598     0:00:23
8   108073540   0:04:47
Yasta
```

which produces all the results you asked for in reasonable times, about two full orders of magnitude (~100x) faster than Emu71, and of course still O(n*log(n)) [ 4:47 / 00:23 = 12.5x ]. A conversion to compiled C# code produces almost two additional orders of magnitude of speed increase.

Quote:

Many thanks to those who contributed solutions, you have some clever and compact solutions, all of which I studied and learned from.

You're welcome and again, thanks for both your interesting challenge and your detailed, well-documented solution.

```
```

Edited: 17 Mar 2012, 12:08 p.m.

Nice chALLENge and even nicer solution!

While I was watching a BBC documentary about Henry the Meek, I decided to try F(100000) on the real HP-42S. The calculator bravely accomplished the task just before the second episode was over.

Let me present an alternative G(N) function for those of you who are in a hurry and can tolerate some error:

```00 { 119-Byte Prgm }
01>LBL "T"
02 STO 00
03 3
04 BASE÷
05 RCL 00
06 2
07 ÷
08 SQRT
09 BASE+
10 RCL 00
11 LN
12 RCL× ST L
13 4.76115738459E-2
14 ×
15 1.29705734564E-1
16 RCL× 00
17 -
18 47.8815
19 -
20 +
21 RCL 00
22 SQRT
23 LASTX
24 LN
25 8.60090977995E-2
26 ×
27 1.50150179296
28 -
29 ×
30 53.456
31 +
32 BASE+
33 .END.
absolute    percent
N         F(N)           G(N)          error       error
1E+02	         42	        43	    +1	   2.381
1E+03	        532	       530	    -2	   0.376
1E+04	       6427	      6426	    -1	   0.016
1E+05	      75243	     75244	    +1	   0.001
1E+06	     861805	    861805	     0	   0.000
1E+07	    9712598	   9712232	  -366	   0.004
1E+08	  108073540	 108074423	  +883	   0.001
1E+09	 1190212172	1190326147     +113975	   0.010
1E+10	12996874312    12999364614    +2490302	   0.019
```

Less than one second on the real HP-42S :-)

Gerson.

I suspect there was some unintended interpretations of the naïve method noted above. Just to be clear when referring to programming and algorithms, "naïve method" referrs to the following definition from http://en.wikipedia.org/wiki/Na%C3%AFve_algorithm#Classification

Quote:
Brute-force or exhaustive search. This is the naïve method method of trying every possible solution to see which is best.

In programming this does not carry any negative connotation- only an acknowledgement that is is a brute-force approach, unoptimized.

Hi, Gerson:

Quote:

Let me present an alternative G(N) function for those of you who are in a hurry and can tolerate some error:

Good idea, here's my take on the subject for the HP-71B:

```1  DESTROY ALL @ OPTION BASE 1 @ DIM C(8) @ READ C @ INPUT K @ FIX 0
2  DATA 2.23467244246E-8,-1.76140742602E-6,5.78605330267E-5,-1.04128697469E-3
3  DATA 1.14617819796E-2,-8.37891596549E-2,2.79672680182,-1.59739297094
4  P=LGT(K) @ N=C(1) @ FOR I=2 TO 8 @ N=N*P+C(I) @ NEXT I @ DISP EXP(N)
```

The results are:

``` P      Exact value  Computed val. Error   % error
----------------------------------------------------
1E02           42           42        0   -0.001
1E03          532          532        0   +0.001
1E04         6427         6427        0   -0.001
1E05        75243        75242       +1   +0.001
1E06       861805       861813       -8   -0.001
1E07      9712598      9712504      +94   +0.001
1E08    108073540    108074589    -1049   -0.001
1E09   1190212172   1190200618   +11554   +0.001
1E10  12996874312  12997000479  -126167   -0.001
```

and of course it runs almost instantaneously.

Best regards from V.

```
```

Hi, Valentin!

Great fit! What about fitting H(N) = F(N) - (N DIV 3 + INT(SQR(N/2))) and then making G(N) = N DIV 3 + INT(SQR(N/2)) + H(N)? Would this result in an even better fit?

Best regards,

Gerson.

Not only a 42S program! It runs with almost no change on my HP-67.

The only thing I needed to change was to replace SIGN on line 8 with 1.

It calculated F(1000) in 2 minutes 26 sec!

**vp

Thank you for this nice challenge.

And felicitations for your brilliant solution. Using Farey series is a much more elegant solution that the crude and brut way to determine coprimes factors from a or b the way I made it!

I am studying your solution with great interest. The use of the coprime factors use from Farey the fractions is really stricky. But I still miss why or how to select range of the serie. Why are FLOOR(SQRT(P)) appropriate ? Is it a kind of magic ? lol
I can't figure out how to relay this with the problem!

I can't resist to present here an old-style user-RPL version using your clever Farey's tricks: It is close-looking to RPL+ version above. Except that I have use more stack manipulations to spare local variables (and speed up runs on my poor old HP-28S).

42 (6s) 532 (20s) 6427 (1min30s)

```« DUP                            @ Preserve P
3 / SQRT FLOOR                 @ Set Order for Recursive Farey Series
1 5 SQRT + 2/                  @ golden ratio
-> P d gr
« 0                            @ Initiate Sum
0 1 1 1                      @ m1 n1 m2 n2 leave and treated in stack
WHILE DUP2 / gr <            @ test  m2/n2 < gr
REPEAT
DUP2 * LAST + SQ SWAP -
P SWAP / FLOOR            @ FLOOR(P/((a+b)²-a*b))
3 PICK d + OVER / FLOOR   @  -> h
3 PICK OVER * 6 ROLL -    @ compute new mt (and remove m1 from stack)
SWAP 3 PICK * 5 ROLL -    @ compute new nt (and remove n1 from stack)
END
4 DROPN                      @ remove m1 n1 m2 n2 from stack (Sum remain on top)
»
»
```

And again, thank to you for this nice challenge and for other forumers to participate. Great Time !

EDIT: remove swapped >> from code listing.

Edited: 19 Mar 2012, 12:00 p.m. after one or more responses were posted

You are most welcome- I enjoyed working on this problem with what time I had.

You raise a good point about how to get the order of the Farey sequence. I had two challenges, first the Sequence normally covers {0..1},but I needed {1..phi} for reasons already discussed.

The second challenge, determining the order = floor(sqrt(p/3)) was not as trivial. In examining the data for the smaller O(N^2) runs, the irreducible coprime factors came up. I noticed that in every case F(N), the boundary closest to unity: k+1/k always stopped at a certain value of k_max (and that no greater denominator existed anywhere in the set):

```  N         k_max
-------------------
100         5
1000       18
10000      57
```

And so on. Since I already came to the same conclusion many others had about the importance of "P/3" in determining the upper limit for side A, It did not take much guess work to identify SQRT(P/3) as an interesting number in the Farey sequence, and the FLOOR function simply ensured the remainder of the denominators were integers.

Quote:
```
« DUP                            @ Preserve P
3 / SQRT FLOOR                 @ Set Order for Recursive Farey Series
1 5 SQRT + 2/                  @ golden ratio
-> P d gr
« 0                            @ Initiate Sum
0 1 1 1                      @ m1 n1 m2 n2 leave and treated in stack
WHILE DUP2 / gr <            @ test  m2/n2 < gr
REPEAT
DUP2 * LAST + SQ SWAP -
P SWAP / FLOOR            @ FLOOR(P/((a+b)²-a*b))
3 PICK d + OVER/ FLOOR    @  -> h
3 PICK OVER * 6 ROLL -    @ compute new mt (and remove m1 from stack)
SWAP 3 PICK * 5 ROLL -    @ compute new nt (and remove n1 from stack)
»
END
4 DROPN                      @ remove m1 n1 m2 n2 from stack (Sum remain on top)
»
```

Are you sure this code will work? Aren't >> and END swapped over?

You are right, the \>> swap here from a bad copy-paste operation!

Corrected listing:

```« DUP                            @ Preserve P
3 / SQRT FLOOR                 @ Set Order for Recursive Farey Series
1 5 SQRT + 2/                  @ golden ratio
-> P d gr
« 0                            @ Initiate Sum
0 1 1 1                      @ m1 n1 m2 n2 leave and treated in stack
WHILE DUP2 / gr <            @ test  m2/n2 < gr
REPEAT
DUP2 * LAST + SQ SWAP -
P SWAP / FLOOR            @ FLOOR(P/((a+b)²-a*b))
3 PICK d + OVER / FLOOR   @  -> h
3 PICK OVER * 6 ROLL -    @ compute new mt (and remove m1 from stack)
SWAP 3 PICK * 5 ROLL -    @ compute new nt (and remove n1 from stack)
END
4 DROPN                      @ remove m1 n1 m2 n2 from stack (Sum remain on top)
»
»
```

Sorry.

Edited: 19 Mar 2012, 12:03 p.m.

```If a <= b <= c
and n/d is the irreducible fraction of b/a,
then the smallest triangle you can construct (subject to b^2 = a*c) is
(a,b,c) = (d^2,n*d,n^2)
because
b = k*n
a = k*d
for k a positive integer, then
b^2 = a*c
k^2 * n^2 = k*d*c
c = k*n^2/d
Since n and d are coprime, k >= d
Since a <= P/3, it follows that
d <= SQRT(P/3)
```

Werner

Werner, A very nice derivation, thank you!

Here is what the code using ISPRIME? and NEXTPRIME instruction may look at.

where ISPRIME? returns true (1) for prime numbers only.
and NEXTPRIME converts its argument to the next prime number (ex. 1 NEXTPRIME is 2, 2 NEXTPRIME is 3, 3 NEXTPRIME is 5, 5 NEXTPRIME is 7 and so on...)

```PROGRAM:                                      @ Comments                    @ Stack   5: 4: 3:  2:  1:
@                             @                            P
« -> P
« 0                                        @ Initiate Sum=0              @                  0 : 1 : P/3
1 P 3 / FOR a                            @  Main Loop                  @                             Sum
2                                    @ Initiate D <- 2
a                                    @          x <- a
1                                    @          i <- 1             @            Sum : D : x :   i
WHILE OVER ISPRIME? NOT              @ loop until x is prime
REPEAT
NEXTPRIME                         @    i' <- nextprime ( i )    @            Sum : D : x :   i
0 3 ROLLD                         @    Initiate  n <- 0         @        Sum : D : n : x :   i
WHILE DUP2 MOD NOT                @    while x MOD i = 0
REPEAT
SWAP OVER / SWAP               @       x <- x/i              @    Sum : D :  n  : x/i :   i
ROT 1 + 3 ROLLD                @       n <- n+1              @    Sum : D : n+1 : x/i :   i
END
4 ROLL OVER 5 ROLL 2 /
CEIL ^ * 3 ROLLD                  @       D = D*i^CEIL(n/2)     @         Sum : D' : x   :   i
END
DROP *                               @   D = D*x                   @                    Sum : D*x
a SWAP /                                                           @                   Sum : a/2D
IF P a 3 5 SQRT + * <
THEN 4 P * a / 3 - SQRT 3 -       @ limited by perimeter        @  Sum : a/2D : SQRT(4P/a-3)-3
ELSE 5 SQRT 1 -                   @ limited by triangle side    @  Sum : a/2D :      SQRT(5)-1
END
* IP 1 + +                           @  Sum' = Sum + 1 + INT(...)  @                         Sum'
NEXT
»
»                                                                           @                         Sum
```

The code is not far from the original one without ISPRIME? and NEXTPRIME.
Time is spare since for large arguments, fewer tests or loops are needed to detemine all cofactor. Especely, no more loop are waste to test for composite (ex. i=9 is always negative because i=3 was alreday tested, but no way !)
Most of the time, D is close to a. With a large a great number of loops are no more waste. And when a is large and prim, no loop at all, immediately D=a is detected.

Time may be lost when ISPRIME? and NEXTPRIME are long time runs. This would be the case with my HP-28S if I have to programm ISPRIME? or NEXTPRIME as USER-RPL. Without built-in instruction using internal speed machine code, no gain.

But still, as demonstrate by Allen, finding cofactor using Farey Series is much more tricky and faster than determining D for each a value, since there is no more loops for cofactors determination or prime testing. The Farey Series directly give coprime b / a pairs and a few scale and end-test is only need :

```« DUP                            @ Preserve P
3 / SQRT FLOOR                 @ Set Order for Recursive Farey Series
1 5 SQRT + 2/                  @ golden ratio
-> P d gr
« 0                            @ Initiate Sum
0 1 1 1                      @ m1 n1 m2 n2 leave and treated in stack
WHILE DUP2 / gr <            @ test  m2/n2 < gr
REPEAT
DUP2 * LAST + SQ SWAP -
P SWAP / FLOOR            @ FLOOR(P/((a+b)²-a*b))
3 PICK d + OVER/ FLOOR    @  -> h
3 PICK OVER * 6 ROLL -    @ compute new mt (and remove m1 from stack)
SWAP 3 PICK * 5 ROLL -    @ compute new nt (and remove n1 from stack)
END
4 DROPN                      @ remove m1 n1 m2 n2 from stack (Sum remain on top)
»
»
```

Edited: 21 Mar 2012, 8:03 a.m.

Thanks for this, C.Ret!

Your points about there being no speed gain if the functions aren't accelerated is understood. I was really just curious if I missed something obvious and the code would simplify significantly.

I tried running it but it gets stuck in an endless loop, with ISPRIME? getting 1 as arguments apparently indefinitely.

I'll try to follow your stack diagram and see where the problem is.

You have to check, 1 ISPRIME? my have return 1 as expected.
You have to check also for 2: 2 ISPRIME? may have return 1 (2 is prime).

Perhaps, an inifinite loop result from NEXTPRIME fonction. I have no idea of how it is working.

My assumption is that :
1 NEXTPRIME returns 2, 2 NEXTPRIME returns 3 etc.

But, depending of impletation, maybe NEXTPRIME returns the closest prime number, resulting in the catastrophic chain reaction:

1 NEXTPRIME returns 1 (as 1 is the closest prime number),

4 NEXTPRIME return 5 (as 5 is the closest prime number - greatest or equal to the argument).

In this case, the modified version:

```PROGRAM:                                      @ Comments                    @ Stack   5: 4: 3:  2:  1:
@                             @                            P
« -> P
« 0                                        @ Initiate Sum=0              @                  0 : 1 : P/3
1 P 3 / FOR a
2                                    @ Initiate D <- 2
a                                    @          x <- a
2                                    @          i <- 2             @            Sum : D : x :   i
WHILE OVER ISPRIME? NOT              @ loop until x is prime
REPEAT
NEXTPRIME                         @    i' <- nextprime ( i )    @            Sum : D : x :   i
0 3 ROLLD                         @    Initiate  n <- 0         @        Sum : D : n : x :   i
WHILE DUP2 MOD NOT                @    while x MOD i = 0
REPEAT
SWAP OVER / SWAP               @       x <- x/i              @    Sum : D :  n  : x/i :   i
ROT 1 + 3 ROLLD                @       n <- n+1              @    Sum : D : n+1 : x/i :   i
END
4 ROLL OVER 5 ROLL 2 /
CEIL ^ * 3 ROLLD                  @       D = D*i^CEIL(n/2)     @         Sum : D' : x   :   i
1 +                               @       inc i                 @
END
DROP *                               @   D = D*x                   @                    Sum : D*x
a SWAP /                                                           @                   Sum : a/2D
IF P a 3 5 SQRT + * <
THEN 4 P * a / 3 - SQRT 3 -       @ limited by perimeter        @  Sum : a/2D : SQRT(4P/a-3)-3
ELSE 5 SQRT 1 -                   @ limited by triangle side    @  Sum : a/2D :      SQRT(5)-1
END
* IP 1 + +                           @  Sum' = Sum + 1 + INT(...)  @                         Sum'
NEXT
»
»                                                                           @                         Sum
```

Hi C.Ret.

Believe it or not, 1 is not a prime... ;-)

You're not alone in thinking that it should be one. Here's what the Wikipedia page about primes has to say about the matter:

Quote:
Most early Greeks did not even consider 1 to be a number,[4] so did not consider it a prime. In the 19th century however, many mathematicians did consider the number 1 a prime. For example, Derrick Norman Lehmer's list of primes up to 10,006,721, reprinted as late as 1956,[5] started with 1 as its first prime.[6] Henri Lebesgue is said to be the last professional mathematician to call 1 prime.[7] Although a large body of mathematical work is also valid when calling 1 a prime, the above fundamental theorem of arithmetic does not hold as stated. For example, the number 15 can be factored as 3 · 5 or 1 · 3 · 5. If 1 were admitted as a prime, these two presentations would be considered different factorizations of 15 into prime numbers, so the statement of that theorem would have to be modified. Furthermore, the prime numbers have several properties that the number 1 lacks, such as the relationship of the number to its corresponding value of Euler's totient function or the sum of divisors function.[8][9]

Anyway, I had suspected that and changed my ISPRIME? function to return true for "1" but that led to stack underrun. I also tried starting the FOR loop with 2.

I shall debug this a little later when I have time.

The 50g (and ND1) behavior for these two functions is:

ISPRIME?: the first int for which 1 is returned is 2

NEXTPRIME: returns 2 for any int smaller than 2; any int is valid input, the next closest prime will be returned

This is obviously hard to develop if you don't actually have the functions built-in.