Pairs with common factors

From Rosetta Code
Pairs with common factors is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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



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


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.

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

Phix

with javascript_semantics
function totient(integer n)
    integer tot = n, i = 2
    while i*i<=n do
        if mod(n,i)=0 then
            while true do
                n /= i
                if mod(n,i)!=0 then exit end if
            end while
            tot -= tot/i
        end if
        i += iff(i=2?1:2)
    end while
    if n>1 then
        tot -= tot/n
    end if
    return tot
end function

constant limit = 1e6
sequence a = repeat(0,limit)
atom sumPhi = 0
for n=1 to limit do
    sumPhi += totient(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

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

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