Pairs with common factors

From Rosetta Code
Task
Pairs with common factors
You are encouraged to solve this task according to the task description, using any language you may know.

Generate the sequence where each term n is the count of the pairs (x,y) with 1 < x < y <= n, that have at least one common prime factor.


For instance, when n = 9, examine the pairs:

   (2,3) (2,4) (2,5) (2,6) (2,7) (2,8) (2,9)
   (3,4) (3,5) (3,6) (3,7) (3,8) (3,9)
   (4,5) (4,6) (4,7) (4,8) (4,9)
   (5,6) (5,7) (5,8) (5,9)
   (6,7) (6,8) (6,9)
   (7,8) (7,9)
   (8,9)

Find all of the pairs that have at least one common prime factor:

   (2,4) (2,6) (2,8) (3,6) (3,9) (4,6) (4,8) (6,8) (6,9)

and count them:

   a(9) = 9


Terms may be found directly using the formula:

   a(n) = n × (n-1) / 2 + 1 - 𝚺{i=1..n} 𝚽(i)

where 𝚽() is Phi; the Euler totient function.


For the term a(p), if p is prime, then a(p) is equal to the previous term.


Task
  • Find and display the first one hundred terms of the sequence.
  • Find and display the one thousandth.


Stretch
  • Find and display the ten thousandth, one hundred thousandth, one millionth.


See also



ALGOL 68

Works with: ALGOL 68G version 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.

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
Output:
    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

ALGOL W

As Algol W is limited to 32 bit integers, shows a(100000) but not a(1000000).

begin % finds integers 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 integers 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  %

    % returns the number of integers k where 1 <= k <= n that are mutually   %
    %                                                             prime to n %
    integer procedure totient ( integer value n ) ;     % algorithm from the %
        if       n < 3 then 1   % 2nd Go Sample in the Totient function task %
        else if  n = 3 then 2
        else begin
            integer t, v, i;
            t := n;
            v := n;
            i := 2;
            while i * i <= v do begin
                if v rem i = 0 then begin
                    while v rem i = 0 do v := v div i;
                    t := t - t div i
                end if_v_ren_i_eq_0 ;
                if i = 2 then i := 1;
                i := i + 2
            end while_ii_le_v ;
            if v > 1 then t - t div v else t
        end totient ;

    integer maxNumber, totientSum, nextToShow;
    maxNumber  := 100000;                 % maximum number of terms required %
    totientSum := 0;                              % sum of the totients 1..n %
    nextToShow := 1000;                   % next power of 10 integer to show %
    for n := 1 until maxNumber do begin
        integer an;
        totientSum := totientSum + totient( n );
        an         := ( ( ( n * ( n - 1 ) ) div 2 ) + 1 ) - totientSum;
        if n <= 100 then begin
            writeon( i_w := 4, s_w := 0, " ", an );
            if n rem 10 = 0 then write()
            end
        else if n = nextToShow then begin
            write( i_w := 1, s_w := 0, "a(", n, "): ", an );
            nextToShow := nextToShow * 10
        end if_n_le_100__n_eq_nextToShow
    end for_n

end.
Output:
    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

C

Translation of: Wren
#include <stdio.h>
#include <stdint.h>
#include <stdbool.h>
#include <stdlib.h>
#include <locale.h>

#define MAX 1000000

bool isPrime(int n) {
    if (n < 2) return false;
    if (n%2 == 0) return n == 2;
    if (n%3 == 0) return n == 3;
    int d = 5;
    while (d*d <= n) {
        if (n%d == 0) return false;
        d += 2;
        if (n%d == 0) return false;
        d += 4;
    }
    return true;
}

uint64_t totient(uint64_t n) {
    uint64_t tot = n, i = 2;
    while (i*i <= n) {
        if (!(n%i)) {
            do {n /= i;} while (!(n%i));
            tot -= tot/i;
        }
        if (i == 2) i = 1;
        i += 2;
    }
    if (n > 1) tot -= tot/n;
    return tot;
}

const char *ord(int c) {
    int m = c % 100;
    if (m >= 4 && m <= 20) return "th";
    m %= 10;
    return (m == 1) ? "st" :
           (m == 2) ? "nd" :
           (m == 3) ? "rd" : "th";
}

int main() {
    uint64_t n, sumPhi = 0, *a = (uint64_t *)calloc(MAX, sizeof(uint64_t));
    int i, limit;
    for (n = 1; n <= MAX; ++n) {
        sumPhi += totient(n);
        if (isPrime(n)) {
            a[n-1] = a[n-2];
        } else {
            a[n-1] = n*(n-1)/2 + 1 - sumPhi;
        }
    }
    setlocale(LC_NUMERIC, "");
    printf("Number of pairs with common factors - first one hundred terms:\n");
    for (i = 0; i < 100; ++i) {
        printf("%'5lu  ", a[i]);
        if (!((i+1)%10)) printf("\n");
    }
    printf("\n");
    for (limit = 1; limit <= MAX; limit *= 10) {
        printf("The %'d%s term: %'lu\n", limit, ord(limit), a[limit-1]);
    }
    free(a);
    return 0;
}
Output:
Number of pairs with common factors - first one hundred terms:
    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  1,006  1,040  1,079  1,095  1,148  1,148  1,195  
1,221  1,262  1,262  1,321  1,341  1,384  1,414  1,461  1,461  1,526  
1,544  1,591  1,623  1,670  1,692  1,755  1,755  1,810  1,848  1,907  

The 1st term: 0
The 10th term: 14
The 100th term: 1,907
The 1,000th term: 195,309
The 10,000th term: 19,597,515
The 100,000th term: 1,960,299,247
The 1,000,000th term: 196,035,947,609

C++

#include <algorithm>
#include <cstdint>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <vector>

std::vector<uint32_t> totients;
std::vector<uint32_t> primes;

void listTotients(const uint32_t& maximum) {
	totients.resize(maximum + 1);
	std::iota(totients.begin(), totients.end(), 0);

	for ( uint32_t i = 2; i <= maximum; ++i ) {
		if ( totients[i] == i ) {
			totients[i] = i - 1;
			for ( uint32_t j = i * 2; j <= maximum; j += i ) {
				totients[j] = ( totients[j] / i ) * ( i - 1 );
			}
		}
	}
}

void listPrimeNumbers(const uint32_t& maximum) {
	const uint32_t halfMaximum = ( maximum + 1 ) / 2;
	std::vector<bool> composite(halfMaximum, false);

	for ( uint32_t i = 1, p = 3; i < halfMaximum; p += 2, ++i ) {
		if ( ! composite[i] ) {
			for ( uint32_t j = i + p; j < halfMaximum; j += p ) {
				composite[j] = true;
			}
		}
	}

	primes.push_back(2);
	for ( uint32_t i = 1, p = 3; i < halfMaximum; p += 2, ++i ) {
		if ( ! composite[i] ) {
			primes.push_back(p);
		}
	}
}

int main() {
	const uint32_t maximum = 1'000'000;
	listPrimeNumbers(maximum);
	listTotients(maximum);
	std::vector<uint64_t> pairsCount(maximum + 1, 0);
	uint64_t totientSum = 0;

	for ( uint64_t number = 1; number <= maximum; ++number ) {
		totientSum += totients[number];
		if ( std::binary_search(primes.begin(), primes.end(), number) ) {
			pairsCount[number] = pairsCount[number - 1];
		} else {
			pairsCount[number] = ( number * ( number - 1 ) >> 1 ) - totientSum + 1;
		}
	}

	std::cout << "The first one hundred terms of the number of pairs with common factors:" << std::endl;
	for ( uint32_t number = 1; number <= 100; ++number ) {
		std::cout << std::setw(4) << pairsCount[number] << ( ( number % 10 == 0 ) ? "\n" : " " );
	}
	std::cout << std::endl;

	uint32_t term = 1;
	while ( term <= maximum ) {
		std::cout << std::left << std::setw(14)
                  << "Term " + std::to_string(term) + ": " << pairsCount[term] << std::endl;
		term *= 10;
	}
}
Output:
The first one hundred terms of the number of pairs with common factors:
   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

Term 1:       0
Term 10:      14
Term 100:     1907
Term 1000:    195309
Term 10000:   19597515
Term 100000:  1960299247
Term 1000000: 196035947609

EasyLang

Translation of: C
func isprim num .
   if i < 2
      return 0
   .
   i = 2
   while i <= sqrt num
      if num mod i = 0
         return 0
      .
      i += 1
   .
   return 1
.
func totient n .
   tot = n
   i = 2
   while i * i <= n
      if n mod i = 0
         repeat
            n /= i
            until n mod i <> 0
         .
         tot -= tot / i
      .
      if i = 2
         i = 1
      .
      i += 2
   .
   if n > 1
      tot -= tot / n
   .
   return tot
.
write "1-100:"
for n = 1 to 1000
   sumPhi += totient n
   if isprim n = 1
      a = ap
   else
      a = n * (n - 1) / 2 + 1 - sumPhi
   .
   if n <= 100
      write " " & a
   .
   ap = a
.
print ""
print "1000: " & a

Factor

Works with: Factor version 0.99 2022-04-03
USING: formatting grouping io kernel math math.functions
math.primes.factors prettyprint ranges sequences
tools.memory.private ;

: totient-sum ( n -- sum )
    [1..b] [ totient ] map-sum ;

: a ( n -- a(n) )
    dup [ 1 - * 2 / ] keep totient-sum - ;

"Number of pairs with common factors - first 100 terms:" print
100 [1..b] [ a commas ] map 10 group simple-table. nl

7 <iota> [ dup 10^ a commas "Term #1e%d: %s\n" printf ] each
Output:
Number of pairs with common factors - first 100 terms:
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   1,006 1,040 1,079 1,095 1,148 1,148 1,195
1,221 1,262 1,262 1,321 1,341 1,384 1,414 1,461 1,461 1,526
1,544 1,591 1,623 1,670 1,692 1,755 1,755 1,810 1,848 1,907

Term #1e0: 0
Term #1e1: 14
Term #1e2: 1,907
Term #1e3: 195,309
Term #1e4: 19,597,515
Term #1e5: 1,960,299,247
Term #1e6: 196,035,947,609

FreeBASIC

Function isPrime(n As Uinteger) As Boolean
    If n Mod 2 = 0 Then Return false
    For i As Uinteger = 3 To Int(Sqr(n))+1 Step 2
        If n Mod i = 0 Then Return false
    Next i
    Return true
End Function

Function Totient(n As Uinteger) As Integer
    Dim As Uinteger tot = n, i = 2
    While i*i <= n
        If n Mod i = 0 Then
            Do
                n /= i
            Loop Until n Mod i <> 0
            tot -= tot/i
        End If
        i += Iif(i = 2, 1, 2)
    Wend
    If n > 1 Then tot -= tot/n
    Return tot
End Function

Dim As Uinteger n, limit = 1e6, sumPhi = 0
Dim As Uinteger a(limit)
For n = 1 To limit
    sumPhi += Totient(n)
    a(n) = Iif(isPrime(n), a(n-1), n * (n - 1) / 2 + 1 - sumPhi)
Next n

Print "Number of pairs with common factors - first one hundred terms:"
Dim As Uinteger j, count = 0
For j = 1 To 100
    count += 1
    Print Using "  ##,###"; a(j);
    If(count Mod 10) = 0 Then Print
Next j

Print !"\nThe 1st term: "; a(1)
Print "The 10th term: "; a(10)
Print "The 100th term: "; a(1e2)
Print "The 1,000th term: "; a(1e3)
Print "The 10,000th term: "; a(1e4)
Print "The 100,000th term: "; a(1e5)
Print "The 1,000,000th term: "; a(1e6)
Sleep
Output:
Number of pairs with common factors - first one hundred terms:
    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   1,006   1,040   1,079   1,095   1,148   1,148   1,195
1,221   1,262   1,262   1,321   1,341   1,384   1,414   1,461   1,461   1,526
1,544   1,591   1,623   1,670   1,692   1,755   1,755   1,810   1,848   1,907

The 1st term: 0
The 10th term: 14
The 100th term: 1907
The 1,000th term: 195309
The 10,000th term: 19597515
The 100,000th term: 1960299247
The 1,000,000th term: 196035947609


Go

Translation of: Wren
Library: Go-rcu
package main

import (
    "fmt"
    "rcu"
)

func totient(n uint64) uint64 {
    tot := n
    i := uint64(2)
    for i*i <= n {
        if n%i == 0 {
            for n%i == 0 {
                n /= i
            }
            tot -= tot / i
        }
        if i == 2 {
            i = 1
        }
        i += 2
    }
    if n > 1 {
        tot -= tot / n
    }
    return tot
}

func ord(c int) string {
    m := c % 100
    if m >= 4 && m <= 20 {
        return "th"
    }
    m %= 10
    switch m {
    case 1:
        return "st"
    case 2:
        return "md"
    case 3:
        return "rd"
    default:
        return "th"
    }
}

func main() {
    const max = 1_000_000
    a := make([]uint64, max)
    sumPhi := uint64(0)
    for n := uint64(1); n <= uint64(max); n++ {
        sumPhi += totient(n)
        if rcu.IsPrime(n) {
            a[n-1] = a[n-2]
        } else {
            a[n-1] = n*(n-1)/2 + 1 - sumPhi
        }
    }
    fmt.Println("Number of pairs with common factors - first one hundred terms:")
    rcu.PrintTable(a[:100], 10, 6, true)
    fmt.Println()
    for limit := 1; limit <= max; limit *= 10 {
        fmt.Printf("The %s%s term: %s\n", rcu.Commatize(limit), ord(limit), rcu.Commatize(a[limit-1]))
    }
}
Output:
Number of pairs with common factors - first one hundred terms:
     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  1,006  1,040  1,079  1,095  1,148  1,148  1,195 
 1,221  1,262  1,262  1,321  1,341  1,384  1,414  1,461  1,461  1,526 
 1,544  1,591  1,623  1,670  1,692  1,755  1,755  1,810  1,848  1,907 

The 1st term: 0
The 10th term: 14
The 100th term: 1,907
The 1,000th term: 195,309
The 10,000th term: 19,597,515
The 100,000th term: 1,960,299,247
The 1,000,000th term: 196,035,947,609

J

For this task, because of the summation of euler totient values, it's more efficient to generate the sequence with a slightly different routine than we would use to compute a single value. Thus:
   (1 _1r2 1r2&p. - +/\@:(5&p:)) 1+i.1e2
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
   (1 _1r2 1r2&p.@{: - +/@:(5&p:)) 1+i.1e3
195309
   (1 _1r2 1r2&p.@{: - +/@:(5&p:)) 1+i.1e4
19597515
   (1 _1r2 1r2&p.@{: - +/@:(5&p:)) 1+i.1e5
1960299247
   (1 _1r2 1r2&p.@{: - +/@:(5&p:)) 1+i.1e6
196035947609

Here, p. calculates a polynomial (1 + (-x)/2 + (x^2)/2 in this example), 5&p: is euler's totient function, @{: modifies the polynomial to only operate on the final element of a sequence, +/ is sum and +/\ is running sum, and 1+i.n is the sequence of numbers 1 through n.

Java

import java.util.ArrayList;
import java.util.Collections;
import java.util.List;

public final class PairsWithCommonFactors {

	public static void main(String[] args) {
		final int maximum = 1_000_000;
		listPrimeNumbers(maximum);
		listTotients(maximum);
		long[] pairsCount = new long[maximum + 1];
		long totientSum = 0;
		
		for ( int number = 1; number <= maximum; number++ ) {
			totientSum += totients[number];
			if ( Collections.binarySearch(primes, number) > 0 ) {
				pairsCount[number] = pairsCount[number - 1];
			} else {
				pairsCount[number] = ( (long) number * ( number - 1 ) >> 1 ) - totientSum + 1;
			}
		}
	
		System.out.println("The first one hundred terms of the number of pairs with common factors:");
		for ( int number = 1; number <= 100; number++ ) {
			System.out.print(String.format("%4d%s", pairsCount[number], ( ( number % 10 == 0 ) ? "\n" : " " )));
		}
		System.out.println();
		
		int term = 1;
		while ( term <= maximum ) {
			System.out.println(String.format("%-14s%s", "Term " + term + ": ", pairsCount[term]));
		    term *= 10;
		}
	}
		
	private static void listTotients(int maximum) {
		totients = new int[maximum + 1];
		for ( int i = 0; i <= maximum; i++ ) { 
			totients[i] = i;
		}
		
		for ( int i = 2; i <= maximum; i++ ) {
			if ( totients[i] == i ) {
				totients[i] = i - 1;
				for ( int j = i * 2; j <= maximum; j += i ) {
					totients[j] = ( totients[j] / i ) * ( i - 1 );
				}
			}
		}
	}

	private static void listPrimeNumbers(int maximum) {
		final int halfMaximum = ( maximum + 1 ) / 2;
		boolean[] composite = new boolean[halfMaximum];
		for ( int i = 1, p = 3; i < halfMaximum; p += 2, i++ ) {
			if ( ! composite[i] ) {
				for ( int j = i + p; j < halfMaximum; j += p ) {
					composite[j] = true;
				}
			}
		}
		
		primes = new ArrayList<Integer>(List.of( 2 ));
		for ( int i = 1, p = 3; i < halfMaximum; p += 2, i++ ) {
			if ( ! composite[i] ) { 
				primes.add(p);
			}
		}
	}
	
	private static int[] totients;
	private static List<Integer> primes;

}
Output:
The first one hundred terms of the number of pairs with common factors:
   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

Term 1:       0
Term 10:      14
Term 100:     1907
Term 1000:    195309
Term 10000:   19597515
Term 100000:  1960299247
Term 1000000: 196035947609

jq

Works with jq and gojq, that is, the C and Go implementations of jq.

If using jq, the definition of `_nwise` can be omitted.

Preliminaries

# For the sake of gojq
def _nwise($n):
  def n: if length <= $n then . else .[0:$n] , (.[$n:] | n) end;
  n;

def lpad($len): tostring | ($len - length) as $l | (" " * $l)[:$l] + .;

def is_prime:
  . as $n
  | if ($n < 2)         then false
    elif ($n % 2 == 0)  then $n == 2
    elif ($n % 3 == 0)  then $n == 3
    elif ($n % 5 == 0)  then $n == 5
    elif ($n % 7 == 0)  then $n == 7
    elif ($n % 11 == 0) then $n == 11
    elif ($n % 13 == 0) then $n == 13
    elif ($n % 17 == 0) then $n == 17
    elif ($n % 19 == 0) then $n == 19
    else
      ($n | sqrt) as $rt
      | 23
      | until( . > $rt or ($n % . == 0); .+2)
      | . > $rt
    end;

# jq optimizes the recursive call of _gcd in the following:
def gcd(a;b):
  def _gcd:
    if .[1] != 0 then [.[1], .[0] % .[1]] | _gcd else .[0] end;
  [a,b] | _gcd ;

def count(s): reduce s as $x (0; .+1);

def totient:
  . as $n
  | count( range(0; .) | select( gcd($n; .) == 1) );

The Task

def sumPhi($max):
  reduce range(1; max+1) as $n ({};
    .sumPhi += ($n|totient)
    | if $n | is_prime
      then .a[$n-1] = .a[$n-2] 
      else
        .a[$n-1] = $n * ($n - 1) / 2 + 1 - .sumPhi
      end ) ;

def limits: [ 1, 10, 1e2, 1e3 ];

"Number of pairs with common factors - first one hundred terms:",
(sumPhi( limits[-1] )
 | (.a[0:100] | _nwise(10) | map(lpad(6)) | join(" ") ),
   ( limits[] as $i
     | "The #\($i) term: \(.a[$i-1])" ) )
Output:
     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
The #1 term: 0
The #10 term: 14
The #100 term: 1907
The #1000 term: 195309

Julia

using Formatting
using Primes

pcf(n) = n * (n - 1) ÷ 2 + 1 - sum(totient, 1:n)

foreach(p -> print(rpad(p[2], 5), p[1] % 20 == 0 ? "\n" : ""), pairs(map(pcf, 1:100)))

for expo in 1:6
    println("The ", format(10^expo, commas = true), "th pair with common factors count is ",
       format(pcf(10^expo), commas = true))
end
Output:
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
The 10th pair with common factors count is 14
The 100th pair with common factors count is 1,907
The 1,000th pair with common factors count is 195,309
The 10,000th pair with common factors count is 19,597,515
The 100,000th pair with common factors count is 1,960,299,247
The 1,000,000th pair with common factors count is 196,035,947,609

Mathematica / Wolfram Language

Translation of: Julia
(* Define the prime counting function (pcf) *)
pcf[n_] := n (n - 1)/2 + 1 - Total[EulerPhi /@ Range[n]]

(* Print pairs of (n, pcf(n)) for n from 1 to 100 *)
Do[
  Module[{pair = pcf[n], pairString},
    pairString = ToString[StringPadRight[ToString[pair], 5]];
    If[Mod[n, 20] == 0,
      WriteString["stdout", pairString <> "\n"],
      WriteString["stdout", pairString, " "]
    ]
  ],
  {n, 1, 100}
]

(* Print the 10^expo-th pair with common factors count *)
Do[
  Print["The ", ToString[NumberForm[10^expo, DigitBlock -> 3]], 
  "th pair with common factors count is ", 
  ToString[NumberForm[pcf[10^expo], DigitBlock -> 3]]],
  {expo, 1, 6}
]
Output:
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 
The 10th pair with common factors count is 14
The 100th pair with common factors count is 1,907
The 1,000th pair with common factors count is 195,309
The 10,000th pair with common factors count is 19,597,515
The 100,000th pair with common factors count is 1,960,299,247
The 1,000,000th pair with common factors count is 196,035,947,609

Maxima

/* Define the prime counting function (pcf) */
pcf(n) := n*(n - 1)/2 + 1 - sum(totient(i), i, 1, n);

/* Print pairs of (n, pcf(n)) for n from 1 to 100 */
for n:1 thru 100 do (
  pcf_n : ev(pcf(n), numer),
  if mod(n, 20) = 0 then (
    printf(true, "~4d~%", pcf_n)
  ) else (
    printf(true, "~4d ", pcf_n)
  )
);

/* Print the 10^expo-th pair with common factors count */
for expo:1 thru 6 do (
  pcf_10expo : ev(pcf(10^expo), numer),
  printf(true, "The ~a-th pair with common factors count is ~a~%", 10^expo, pcf_10expo)
);
Output:
   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
The 10-th pair with common factors count is 14
The 100-th pair with common factors count is 1907
The 1000-th pair with common factors count is 195309
The 10000-th pair with common factors count is 19597515
The 100000-th pair with common factors count is 1960299247
The 1000000-th pair with common factors count is 196035947609


Nim

Translation of: Wren
import std/[sequtils, strutils]

func isPrime(n: Natural): bool =
  if n < 2: return false
  if (n and 1) == 0: return n == 2
  if n mod 3 == 0: return n == 3
  var k = 5
  while k * k <= n:
    if n mod k == 0 or n mod (k + 2) == 0: return false
    inc k, 6
  result = true

func totient(n: Natural): int =
  var n = n
  result = n
  var i = 2
  while i * i <= n:
    if n mod i == 0:
      while n mod i == 0:
        n = n div i
      dec result, result div i
    if i == 2: i = 1
    inc i, 2
  if n > 1:
    dec result, result div n

const Max = 1_000_000
var a: array[Max, int]
var sumPhi = 0
for n in 1..Max:
  inc sumPhi, n.totient
  if n.isPrime:
    a[n - 1] = a[n - 2]
  else:
    a[n-1] = n * (n - 1) shr 1 + 1 - sumPhi

echo "Number of pairs with common factors - First one hundred terms:"
for n in countup(0, 99, 10):
  echo a[n..n+9].mapIt(align($it, 4)).join(" ")
echo()
var limit = 1
while limit <= Max:
  echo "The $1th term: $2".format(insertSep($limit), insertSep($a[limit-1]))
  limit *= 10
Output:
Number of pairs with common factors - First one hundred terms:
   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

The 1th term: 0
The 10th term: 14
The 100th term: 1_907
The 1_000th term: 195_309
The 10_000th term: 19_597_515
The 100_000th term: 1_960_299_247
The 1_000_000th term: 196_035_947_609

Pascal

Free Pascal

modified Perfect_totient_numbers

program PairsWithCommonFactors;
{$IFdef FPC} {$MODE DELPHI} {$Optimization ON,ALL}{$IFEND}
{$IFdef Windows} {$APPTYPE CONSOLE}{$IFEND}
const
  cLimit = 1000*1000*1000;
//global
type
  TElem= Uint64;
  tpElem = pUint64;

  myString = String[31];

var
  TotientList : array of TElem;
  Sieve : Array of byte;
  
function Numb2USA(n:Uint64):myString;
const
//extend s by the count of comma to be inserted
  deltaLength : array[0..24] of byte =
    (0,0,0,0,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7,7,7);
var
  pI :pChar;
  i,j : NativeInt;
Begin
  str(n,result);
  i := length(result);
 //extend s by the count of comma to be inserted
// j := i+ (i-1) div 3;
  j := i+deltaLength[i];
  if i<> j then
  Begin
    setlength(result,j);
    pI := @result[1];
    dec(pI);
    while i > 3 do
    Begin
       //copy 3 digits
       pI[j] := pI[i];
       pI[j-1] := pI[i-1];
       pI[j-2] := pI[i-2];
       // insert comma
       pI[j-3] := ',';
       dec(i,3);
       dec(j,4);
    end;
  end;
end;

procedure SieveInit(svLimit:NativeUint);
var
  pSieve:pByte;
  i,j,pr :NativeUint;
Begin
  svlimit := (svLimit+1) DIV 2;
  setlength(sieve,svlimit+1);
  pSieve := @Sieve[0];
  For i := 1 to svlimit do
  Begin
    IF pSieve[i]= 0 then
    Begin
      pr := 2*i+1;
      j := (sqr(pr)-1) DIV 2;
      IF  j> svlimit then
        BREAK;
      repeat
        pSieve[j]:= 1;
        inc(j,pr);
      until j> svlimit;
    end;
  end;
  pr := 0;
  j := 0;
  For i := 1 to svlimit do
  Begin
    IF pSieve[i]= 0 then
    Begin
      pSieve[j] := i-pr;
      inc(j);
      pr := i;
    end;
  end;
  setlength(sieve,j);
end;

procedure TotientInit(len: NativeUint);
var
  pTotLst : tpElem;
  pSieve  : pByte;
  i: NativeInt;
  p,j,k,svLimit : NativeUint;
Begin
  SieveInit(len);
  setlength(TotientList,len+12);
  pTotLst := @TotientList[0];

//Fill totient with simple start values for odd and even numbers
//and multiples of 3
  j := 1;
  k := 1;// k == j DIV 2
  p := 1;// p == j div 3;
  repeat
    pTotLst[j] := j;//1
    pTotLst[j+1] := k;//2 j DIV 2; //2
    inc(k);
    inc(j,2);
    pTotLst[j] := j-p;//3
    inc(p);
    pTotLst[j+1] := k;//4  j div 2
    inc(k);
    inc(j,2);
    pTotLst[j] := j;//5
    pTotLst[j+1] := p;//6   j DIV 3 <=  (div 2) * 2 DIV/3
    inc(j,2);
    inc(p);
    inc(k);
  until j>len+6;

//correct values of totient by prime factors
  svLimit := High(sieve);
  p := 3;// starting after 3
  pSieve := @Sieve[svLimit+1];
  i := -svlimit;
  repeat
    p := p+2*pSieve[i];
    j := p;
    while j <= cLimit do
    Begin
      k:= pTotLst[j];
      pTotLst[j]:= k-(k DIV p);
      inc(j,p);
    end;
    inc(i);
  until i=0;
  //primes not needed anymore
  setlength(sieve,0);
end;

procedure CountOfPairs(len: NativeUint);
var
  pTotLst : tpElem;
  i,a_n,sum,Totient: tElem;
Begin
  pTotLst := @TotientList[0];
  sum := 1;
  a_n := 2; // sums to i*(i-1)/2 +1
  For i := 2 to len do
  Begin
    Totient := pTotLst[i];// relict for print data
    sum += Totient;
    pTotLst[i] := a_n-sum;
    a_n += i;
  end;
  TotientList[1] := 0;
end;

var
  i,k : NativeUint;
Begin
  TotientInit(climit);
  CountOfPairs(climit);
  i := 1;
  Repeat
    For k := 9 downto 0 do
    begin
      write(TotientList[i]:6);
      inc(i);
    end;
    writeln;
  until i>99;
  writeln;
  writeln('Some values #');
  i := 10;
  repeat
    writeln(Numb2USA(i):13,Numb2USA(TotientList[i]):25);
    i *= 10;
  until i > climit;
end.
Output:
     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

Some values #
           10                       14
          100                    1,907
        1,000                  195,309
       10,000               19,597,515
      100,000            1,960,299,247
    1,000,000          196,035,947,609
   10,000,000       19,603,638,572,759
  100,000,000    1,960,364,433,634,093
1,000,000,000  196,036,448,326,991,587
real    0m23,577s

Perl

Library: ntheory
use v5.36;
use ntheory 'factor';
use List::Util qw<sum product uniq max>;

sub comma { reverse ((reverse shift) =~ s/(.{3})/$1,/gr) =~ s/^,//r }
sub table ($c, @V) { my $t = $c * (my $w = 2 + max map {length} @V); ( sprintf( ('%'.$w.'s')x@V, @V) ) =~ s/.{1,$t}\K/\n/gr }

my($max, @phi, @n_pairs) = (100, 0);
for my $t (1..$max) { push @phi, $t * product map { 1 - 1/$_ } uniq factor($t) }
push @n_pairs, comma $_ * ($_ - 1) / 2 + 1 - sum @phi[1..$_] for 1..$max;

say 'Number of pairs with common factors - first one hundred terms:';
say table 10, @n_pairs;
Output:
Number of pairs with common factors - first one hundred terms:
      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  1,006  1,040  1,079  1,095  1,148  1,148  1,195
  1,221  1,262  1,262  1,321  1,341  1,384  1,414  1,461  1,461  1,526
  1,544  1,591  1,623  1,670  1,692  1,755  1,755  1,810  1,848  1,907

Phix

with javascript_semantics
constant limit = 1e6
sequence a = repeat(0,limit)
atom sumPhi = 0
for n=1 to limit do
    sumPhi += phi(n)
    if is_prime(n) then
        a[n] = a[n-1]
    else
        a[n] = n * (n - 1) / 2 + 1 - sumPhi
    end if
end for

string j = join_by(a[1..100],1,10,fmt:="%,5d") 
printf(1,"Number of pairs with common factors - first one hundred terms:\n%s\n",j)
for l in {1, 10, 1e2, 1e3, 1e4, 1e5, 1e6} do
    printf(1,"%22s term: %,d\n", {proper(ordinal(l),"SENTENCE"), a[l]})
end for
Output:
Number of pairs with common factors - first one hundred terms:
    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   1,006   1,040   1,079   1,095   1,148   1,148   1,195
1,221   1,262   1,262   1,321   1,341   1,384   1,414   1,461   1,461   1,526
1,544   1,591   1,623   1,670   1,692   1,755   1,755   1,810   1,848   1,907

                 First term: 0
                 Tenth term: 14
         One hundredth term: 1,907
        One thousandth term: 195,309
        Ten thousandth term: 19,597,515
One hundred thousandth term: 1,960,299,247
         One millionth term: 196,035,947,609

Quackery

totient is defined at Totient function#Quackery.

  [] 0
  1000 times
    [ i^ 1+ totient +
      i^ 1+ dup 1 - * 2 / 1+
      over -
      swap dip join ]
  drop
  dup -1 peek swap
  100 split drop
  say "First 100 terms:"
  [] swap witheach
     [ number$ nested join ]
  60 wrap$
  cr cr say "1000th term: " echo
Output:
First 100 terms:
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

1000th term: 195309

Raku

use Prime::Factor;
use Lingua::EN::Numbers;

my \𝜑 = 0, |(1..*).hyper.map: -> \t { t × [×] t.&prime-factors.unique.map: { 1 - 1/$_ } }

sub pair-count (\n) {  n × (n - 1) / 2 + 1 - sum 𝜑[1..n] }

say "Number of pairs with common factors - first one hundred terms:\n"
    ~ (1..100).map(&pair-count).batch(10)».&comma».fmt("%6s").join("\n") ~ "\n";

for ^7 { say (my $i = 10 ** $_).&ordinal.tc.fmt("%22s term: ") ~ $i.&pair-count.&comma }
Output:
Number of pairs with common factors - first one hundred terms:
     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  1,006  1,040  1,079  1,095  1,148  1,148  1,195
 1,221  1,262  1,262  1,321  1,341  1,384  1,414  1,461  1,461  1,526
 1,544  1,591  1,623  1,670  1,692  1,755  1,755  1,810  1,848  1,907

                 First term: 0
                 Tenth term: 14
         One hundredth term: 1,907
        One thousandth term: 195,309
        Ten thousandth term: 19,597,515
One hundred thousandth term: 1,960,299,247
         One millionth term: 196,035,947,609

RPL

The PHI function is defined at Totient function To save time, Σφ(j) are stored in a global variable named ΣPHI, which must be initialized at { 1 } once.

Translation of: C
Works with: Halcyon Calc version 4.2.8
RPL code Comment
∑PHI SIZE 
  IF DUP2 ≤ THEN DROP 
  ELSE 
     '∑PHI' OVER GET SWAP 1 + ∑PHI SWAP 4 ROLL FOR j
        SWAP j PHI + SWAP OVER + NEXT
     '∑PHI' STO DROP ∑PHI SIZE END
 '∑PHI' SWAP GET 
≫ '→∑PHI' STO

≪ DUP DUP 1 - * 2 / 1 + SWAP →∑PHI - ≫ 'A185670' STO
'→∑PHI' ( n -- φ(n) )
if n... 
...is not already in ∑PHI
  initialize stack and loop, for size(∑PHI)+1 to n
  sum += φ(n), append sum to list 
store list into ∑PHI
read ∑PHI[n]


'A185670' ( n -- n*(n-1)/2 + 1 - Σφ(j) )
≪ { } 1 100 FOR j j A185670 + NEXT ≫ EVAL
1000 A185670
Works with: HP version 49
« DUP DUP 1 - 2 * / 1 +
  'j' 4 ROLL 'EULER(j)' ∑ - R→I
» 'A185670' STO

« « n A185670 » 1 100 1 SEQ
  1000 A185670
» 'TASK' STO
Output:
2: { 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 }
1: 195309

Ruby

require "prime"

 def 𝜑(n) = n.prime_division.inject(1) {|res, (pr, exp)| res *= (pr-1) * pr**(exp-1) } 
 def a(n) = n*(n-1)/2 + 1 - (1..n).sum{|i| 𝜑(i)}

 puts "Number of pairs with common factors - first 100 terms: " 
 puts (1..100).map{|n| a(n) }.join(", ")
 (1..6).each{|e| puts "Term #1e#{e}: #{a(10**e)}"}
Output:
Number of pairs with common factors - first 100 terms: 
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
1e1: 14
1e2: 1907
1e3: 195309
1e4: 19597515
1e5: 1960299247
1e6: 196035947609

Wren

Library: Wren-math
Library: Wren-fmt
import "./math" for Int
import "./fmt" for Fmt

var totient = Fn.new { |n|
    var tot = n
    var i = 2
    while (i*i <= n) {
        if (n%i == 0) {
            while(n%i == 0) n = (n/i).floor
            tot = tot - (tot/i).floor
        }
        if (i == 2) i = 1
        i = i + 2
    }
    if (n > 1) tot = tot - (tot/n).floor
    return tot
}

var max = 1e6
var a = List.filled(max, 0)
var sumPhi = 0
for (n in 1..max) {
    sumPhi = sumPhi + totient.call(n)
    if (Int.isPrime(n)) {
        a[n-1] = a[n-2]
    } else {
        a[n-1] = n * (n - 1) / 2 + 1 - sumPhi
    }
}

System.print("Number of pairs with common factors - first one hundred terms:")
Fmt.tprint("$,5d ", a[0..99], 10)
System.print()
var limits = [1, 10, 1e2, 1e3, 1e4, 1e5, 1e6]
for (limit in limits) {
    Fmt.print("The $,r term: $,d", limit, a[limit-1])
}
Output:
Number of pairs with common factors - first one hundred terms:
    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  1,006  1,040  1,079  1,095  1,148  1,148  1,195  
1,221  1,262  1,262  1,321  1,341  1,384  1,414  1,461  1,461  1,526  
1,544  1,591  1,623  1,670  1,692  1,755  1,755  1,810  1,848  1,907  

The 1st term: 0
The 10th term: 14
The 100th term: 1,907
The 1,000th term: 195,309
The 10,000th term: 19,597,515
The 100,000th term: 1,960,299,247
The 1,000,000th term: 196,035,947,609