Pairs with common factors: Difference between revisions

=={{header|ALGOL 68}}==

{{works with|ALGOL 68G|3 - tested with release 3.0.3.win32}}

Should work with other implementations of Algol 68 where LONG INT is at least 64 bits.

Unfortunately, Algol 68G version 2 runs out of memory sometime after the 100 000th element, however version 3 has no such problem.

<syntaxhighlight lang="algol68">

BEGIN # finds elements of the sequence a(n) where a(n) is number of pairs #

# (x,y) where 1 < x < y <= n that have at least one common prime #

# factor. The sequence elements can be calculated by: #

# a(n) = n(n-1)/2 + 1 - sum i = 1..n of phi(i) where phi is Euler's #

# totient function #

MODE ELEMENT = LONG INT; # integral type large enough for a(1 000 000) #

# returns the number of integers k where 1 <= k <= n that are mutually #

# prime to n #

PROC totient = ( ELEMENT n )ELEMENT: # algorithm from the second #

IF n < 3 THEN 1 # Go Sample in the Totient function task #

ELIF n = 3 THEN 2

ELSE

ELEMENT result := n;

ELEMENT v := n;

ELEMENT i := 2;

WHILE i * i <= v DO

IF v MOD i = 0 THEN

WHILE v MOD i = 0 DO v OVERAB i OD;

result -:= result OVER i

FI;

IF i = 2 THEN

i := 1

FI;

i +:= 2

OD;

IF v > 1 THEN result -:= result OVER v FI;

result

FI # totient # ;

INT max number = 1 000 000; # maximum number of terms required #

ELEMENT totient sum := 0; # sum of the totients 1..n #

INT next to show := 1 000; # next power of 10 element to show #

ELEMENT n := 0;

TO max number DO

n +:= 1;

totient sum +:= totient( n );

ELEMENT an = ( ( ( n * ( n - 1 ) ) OVER 2 ) + 1 ) - totient sum;

IF n <= 100 THEN

print( ( " ", whole( an, -4 ) ) );

IF n MOD 10 = 0 THEN print( ( newline ) ) FI

ELIF n = next to show THEN

print( ( "a(", whole( n, 0 ), "): ", whole( an, 0 ), newline ) );

next to show *:= 10

FI

OD

END

</syntaxhighlight>

{{out}}

<pre>

0 0 0 1 1 4 4 7 9 14

14 21 21 28 34 41 41 52 52 63

71 82 82 97 101 114 122 137 137 158

158 173 185 202 212 235 235 254 268 291

291 320 320 343 363 386 386 417 423 452

470 497 497 532 546 577 597 626 626 669

669 700 726 757 773 818 818 853 877 922

922 969 969 1006 1040 1079 1095 1148 1148 1195

1221 1262 1262 1321 1341 1384 1414 1461 1461 1526

1544 1591 1623 1670 1692 1755 1755 1810 1848 1907

a(1000): 195309

a(10000): 19597515

a(100000): 1960299247

a(1000000): 196035947609

</pre>