Greatest prime dividing the n-th cubefree number: Difference between revisions

From Rosetta Code
Content added Content deleted
(Added Python)
(Added Dart)
Line 24: Line 24:
;Reference
;Reference
* [https://oeis.org/A370833 OEIS sequence: A370833: a(n) is the greatest prime dividing the n-th cubefree number, for n >= 2; a(1)=1.]
* [https://oeis.org/A370833 OEIS sequence: A370833: a(n) is the greatest prime dividing the n-th cubefree number, for n >= 2; a(1)=1.]

=={{header|Dart}}==
{{trans|Wren}}
<syntaxhighlight lang="dart">import 'dart:math';

void main() {
List<int> res = [1];
int count = 1;
int i = 2;
int lim1 = 100;
int lim2 = 1000;
double max = 1e7;
var t0 = DateTime.now();

while (count < max) {
bool cubeFree = false;
List<int> factors = primeFactors(i);
if (factors.length < 3) {
cubeFree = true;
} else {
cubeFree = true;
for (int i = 2; i < factors.length; i++) {
if (factors[i - 2] == factors[i - 1] && factors[i - 1] == factors[i]) {
cubeFree = false;
break;
}
}
}
if (cubeFree) {
if (count < lim1) res.add(factors.last);
count += 1;
if (count == lim1) {
print("First $lim1 terms of a[n]:");
print(res.take(lim1).join(', '));
print("");
} else if (count == lim2) {
print("The $count term of a[n] is ${factors.last}");
lim2 *= 10;
}
}
i += 1;
}
print("${DateTime.now().difference(t0).inSeconds} sec.");
}

List<int> primeFactors(int n) {
List<int> factors = [];
while (n % 2 == 0) {
factors.add(2);
n ~/= 2;
}
for (int i = 3; i <= sqrt(n); i += 2) {
while (n % i == 0) {
factors.add(i);
n ~/= i;
}
}
if (n > 2) {
factors.add(n);
}
return factors;
}</syntaxhighlight>
{{out}}
<pre>First 100 terms of a[n]:
1, 2, 3, 2, 5, 3, 7, 3, 5, 11, 3, 13, 7, 5, 17, 3, 19, 5, 7, 11, 23, 5, 13, 7, 29, 5, 31, 11, 17, 7, 3, 37, 19, 13, 41, 7, 43, 11, 5, 23, 47, 7, 5, 17, 13, 53, 11, 19, 29, 59, 5, 61, 31, 7, 13, 11, 67, 17, 23, 7, 71, 73, 37, 5, 19, 11, 13, 79, 41, 83, 7, 17, 43, 29, 89, 5, 13, 23, 31, 47, 19, 97, 7, 11, 5, 101, 17, 103, 7, 53, 107, 109, 11, 37, 113, 19, 23, 29, 13, 59

The 1000 term of a[n] is 109
The 10000 term of a[n] is 101
The 100000 term of a[n] is 1693
The 1000000 term of a[n] is 1202057
The 10000000 term of a[n] is 1202057</pre>


=={{header|FreeBASIC}}==
=={{header|FreeBASIC}}==

Revision as of 01:02, 6 March 2024

Greatest prime dividing the n-th cubefree number 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.
Definitions

A cubefree number is a positive integer whose prime factorization does not contain any third (or higher) power factors. If follows that all primes are trivially cubefree and the first cubefree number is 1 because it has no prime factors.

Let a[n] be the greatest prime dividing the n-th cubefree number for n >= 2. By convention, let a[1] = 1 even though the first cubefree number, 1, has no prime factors.


Examples

a[2] is clearly 2 because it is the second cubefree number and also prime. The fourth cubefree number is 4 and it's highest prime factor is 2, hence a[4] = 2.


Task

Compute and show on this page the first 100 terms of a[n]. Also compute and show the 1,000th, 10,000th and 100,000th members of the sequence.


Stretch

Compute and show the 1 millionth and 10 millionth terms of the sequence.

This may take a while for interpreted languages.


Reference

Dart

Translation of: Wren
import 'dart:math';

void main() {
  List<int> res = [1];
  int count = 1;
  int i = 2;
  int lim1 = 100;
  int lim2 = 1000;
  double max = 1e7;
  var t0 = DateTime.now();

  while (count < max) {
    bool cubeFree = false;
    List<int> factors = primeFactors(i);
    if (factors.length < 3) {
      cubeFree = true;
    } else {
      cubeFree = true;
      for (int i = 2; i < factors.length; i++) {
        if (factors[i - 2] == factors[i - 1] && factors[i - 1] == factors[i]) {
          cubeFree = false;
          break;
        }
      }
    }
    if (cubeFree) {
      if (count < lim1) res.add(factors.last);
      count += 1;
      if (count == lim1) {
        print("First $lim1 terms of a[n]:");
        print(res.take(lim1).join(', '));
        print("");
      } else if (count == lim2) {
        print("The $count term of a[n] is ${factors.last}");
        lim2 *= 10;
      }
    }
    i += 1;
  }
  print("${DateTime.now().difference(t0).inSeconds} sec.");
}

List<int> primeFactors(int n) {
  List<int> factors = [];
  while (n % 2 == 0) {
    factors.add(2);
    n ~/= 2;
  }
  for (int i = 3; i <= sqrt(n); i += 2) {
    while (n % i == 0) {
      factors.add(i);
      n ~/= i;
    }
  }
  if (n > 2) {
    factors.add(n);
  }
  return factors;
}
Output:
First 100 terms of a[n]:
1, 2, 3, 2, 5, 3, 7, 3, 5, 11, 3, 13, 7, 5, 17, 3, 19, 5, 7, 11, 23, 5, 13, 7, 29, 5, 31, 11, 17, 7, 3, 37, 19, 13, 41, 7, 43, 11, 5, 23, 47, 7, 5, 17, 13, 53, 11, 19, 29, 59, 5, 61, 31, 7, 13, 11, 67, 17, 23, 7, 71, 73, 37, 5, 19, 11, 13, 79, 41, 83, 7, 17, 43, 29, 89, 5, 13, 23, 31, 47, 19, 97, 7, 11, 5, 101, 17, 103, 7, 53, 107, 109, 11, 37, 113, 19, 23, 29, 13, 59

The 1000 term of a[n] is 109
The 10000 term of a[n] is 101
The 100000 term of a[n] is 1693
The 1000000 term of a[n] is 1202057
The 10000000 term of a[n] is 1202057

FreeBASIC

Translation of: Wren

Without using external libraries, it takes about 68 seconds to run on my system (Core i5). 

Dim Shared As Integer factors()
Dim As Integer res(101)
res(0) = 1
Dim As Integer count = 1
Dim As Integer j, i = 2
Dim As Integer lim1 = 100
Dim As Integer lim2 = 1000
Dim As Integer max = 1e7
Dim As Integer cubeFree = 0

Sub primeFactors(n As Integer, factors() As Integer)
    Dim As Integer i = 2, cont = 2
    While (i * i <= n)
        While (n Mod i = 0)
            n /= i
            cont += 1
            Redim Preserve factors(1 To cont)
            factors(cont) = i
        Wend
        i += 1
    Wend
    If (n > 1) Then
        cont += 1
        Redim Preserve factors(1 To cont)
        factors(cont) = n
    End If
End Sub

While (count < max)
    primeFactors(i, factors())
    If (Ubound(factors) < 3) Then
        cubeFree = 1
    Else
        cubeFree = 1
        For j = 2 To Ubound(factors)
            If (factors(j-2) = factors(j-1) And factors(j-1) = factors(j)) Then
                cubeFree = 0
                Exit For
            End If
        Next j
    End If
    
    If (cubeFree = 1) Then
        If (count < lim1) Then
            res(count) = factors(Ubound(factors))
        End If
        count += 1
        If (count = lim1) Then
            Print "First "; lim1; " terms of a[n]:"
            For j = 1 To lim1
                Print Using "####"; res(j-1);
                If (j Mod 10 = 0) Then Print
            Next j
            Print
        Elseif (count = lim2) Then
            Print "The"; count; " term of a[n] is"; factors(Ubound(factors))
            lim2 *= 10
        End If
    End If
    i += 1
Wend

Sleep
Output:
First 100 terms of a[n]:
  1   2   3   2   5   3   7   3   5  11
  3  13   7   5  17   3  19   5   7  11
 23   5  13   7  29   5  31  11  17   7
  3  37  19  13  41   7  43  11   5  23
 47   7   5  17  13  53  11  19  29  59
  5  61  31   7  13  11  67  17  23   7
 71  73  37   5  19  11  13  79  41  83
  7  17  43  29  89   5  13  23  31  47
 19  97   7  11   5 101  17 103   7  53
107 109  11  37 113  19  23  29  13  59

The 1000th term of a[n] is 109
The 10000th term of a[n] is 101
The 100000th term of a[n] is 1693
The 1000000th term of a[n] is 1202057
The 10000000th term of a[n] is 1202057

Phix

Quite possibly flawed, but it does finish in the blink of an eye. 

with javascript_semantics

-- Note this routine had a major hiccup at 27000 = 2^3*3^3*5^3 and another was predicted at 9261000 = 27000*7^3.
-- The "insufficient kludge" noted below makes it match Wren over than limit, but quite possibly only by chance.
function cubes_before(atom n)
    -- nb: if n is /not/ cube-free it /is/ included in the result.
    -- nth := n-cubes_before(n) means n is the nth cube-free integer,
    --                but only if cubes_before(n-1)==cubes_before(n),
    -- otherwise cubes_before(cubicate) isn't really very meaningful.
    integer r = 0
    bool xpm = true -- extend prior multiples
    sequence pm = {}
    for i=1 to n do -- (or eg floor(cbrt(n)))
        atom p3 = power(get_prime(i),3)
        if p3>n then exit end if
        integer k = floor(n/p3)
        for mask=1 to power(2,length(pm))-1 do
            integer m = mask, mi = 0
            atom kpm = p3  
            while m do  
                mi += 1
                if odd(m) then
                    kpm *= pm[mi]
                end if
                m = floor(m/2)
            end while
            if kpm>n then
                if count_bits(mask)=1 then
                    xpm = false
                    pm = pm[1..mi-1]
                    exit
                end if
            else
                -- subtract? already counted multiples..
                integer l = floor(n/kpm)
-- DEV insufficient kludge... (probably)
--if count_bits(mask)=1 then
if odd(count_bits(mask)) then -- match Wren
                k -= l
else
                k += l
end if
            end if
        end for
        r += k
        if xpm then
            pm &= p3
        end if
    end for
    return r
end function

function cube_free(atom nth)
    -- get the nth cube-free integer
    atom lo = nth, hi = lo*2, mid, cb, k
    while hi-cubes_before(hi)<nth do
        lo = hi
        hi = lo*2
    end while
    -- bisect until we have a possible...
    while true do
        mid = floor((lo+hi)/2)
        cb = cubes_before(mid)
        k = mid-cb
        if k=nth then
            -- skip omissions
            while cubes_before(mid-1)!=cb do
                mid -= 1
                cb -= 1
            end while
            exit
        elsif k<nth then
            lo = mid
        else
            hi = mid
        end if
    end while 
    return mid
end function

function A370833(atom n)
    if n=1 then return 1 end if
    atom nth = cube_free(n)
    sequence f = prime_powers(nth)
    return f[$][1]
end function

atom t0 = time()
printf(1,"First 100 terms of a[n]:\n%s\n",join_by(apply(tagset(100),A370833),1,10," ",fmt:="%3d"))
for n in sq_power(10,tagset(7,3)) do
    printf(1,"The %,dth term of a[n] is %,d.\n",{n,A370833(n)})
end for
?elapsed(time()-t0)
Output:
First 100 terms of a[n]:
  1   2   3   2   5   3   7   3   5  11
  3  13   7   5  17   3  19   5   7  11
 23   5  13   7  29   5  31  11  17   7
  3  37  19  13  41   7  43  11   5  23
 47   7   5  17  13  53  11  19  29  59
  5  61  31   7  13  11  67  17  23   7
 71  73  37   5  19  11  13  79  41  83
  7  17  43  29  89   5  13  23  31  47
 19  97   7  11   5 101  17 103   7  53
107 109  11  37 113  19  23  29  13  59

The 1,000th term of a[n] is 109.
The 10,000th term of a[n] is 101.
The 100,000th term of a[n] is 1,693.
The 1,000,000th term of a[n] is 1,202,057.
The 10,000,000th term of a[n] is 1,202,057.
"0.1s"

Python

Library: SymPy
Works with: Python version 3.x
Translation of: FreeBASIC

This takes under 12 minutes to run on my system (Core i5). 

#!/usr/bin/python

from sympy import primefactors

res = [1]
count = 1
i = 2
lim1 = 100
lim2 = 1000
max = 1e7

while count < max:
    cubeFree = False
    factors = primefactors(i)
    if len(factors) < 3:
        cubeFree = True
    else:
        cubeFree = True
        for j in range(2, len(factors)):
            if factors[j-2] == factors[j-1] and factors[j-1] == factors[j]:
                cubeFree = False
                break
    if cubeFree:
        if count < lim1:
            res.append(factors[-1])
        count += 1
        if count == lim1:
            print("First {} terms of a[n]:".format(lim1))
            for k in range(0, len(res), 10):
                print(" ".join(map(str, res[k:k+10])))
            print("")
        elif count == lim2:
            print("The {} term of a[n] is {}".format(count, factors[-1]))
            lim2 *= 10
    i += 1
Output:
Similar to FreeBASIC entry.

RPL

I confirm it does take a while for an interpreted language like RPL. Getting the 100,000th term is likely to be a question of hours, even on an emulator.

Works with: HP version 49
« 1 { } DUP2 + → n res2 res1
  « 2
    WHILE 'n' INCR 10000 ≤ REPEAT
      WHILE DUP FACTORS DUP 1 « 3 < NSUB 2 MOD NOT OR » DOSUBS ΠLIST NOT 
      REPEAT DROP 1 + END
      HEAD
      CASE 
         n 100 ≤      THEN 'res1' OVER STO+ END
         n LOG FP NOT THEN 'res2' OVER STO+ END
      END
      DROP 1 +
    END
    DROP res1 res2
» » 'TASK' STO
Output:
2: { 1 2 3 2 5 3 7 3 5 11 3 13 7 5 17 3 19 5 7 11 23 5 13 7 29 5 31 11 17 7 3 37 19 13 41 7 43 11 5 23 47 7 5 17 13 53 11 19 29 59 5 61 31 7 13 11 67 17 23 7 71 73 37 5 19 11 13 79 41 83 7 17 43 29 89 5 13 23 31 47 19 97 7 11 5 101 17 103 7 53 107 109 11 37 113 19 23 29 13 59 }
1: { 109 101 }

Wren

Library: Wren-math
Library: Wren-fmt

This takes just under 4 minutes to run on my system (Core i7).

Although it looks odd that the 1 millionth and 10 millionth terms are the same (and I have no independent source to check against), it would appear that the 1 millionth cubefree number is 1,202,057 which is prime and the 10 millionth such number is coincidentally 10 times this. 

import "./math" for Int
import "./fmt" for Fmt

var res = [1]
var count = 1
var i = 2
var lim1 = 100
var lim2 = 1000
var max = 1e7
while (count < max) {
    var cubeFree = false
    var factors = Int.primeFactors(i)
    if (factors.count < 3) {
        cubeFree = true
    } else {
        cubeFree = true
        for (i in 2...factors.count) {
            if (factors[i-2] == factors[i-1] && factors[i-1] == factors[i]) {
                cubeFree = false
                break
            }
        }
    }
    if (cubeFree) {
        if (count < lim1) res.add(factors[-1])
        count = count + 1
        if (count == lim1) {
            System.print("First %(lim1) terms of a[n]:")
            Fmt.tprint("$3d", res, 10)
        } else if (count == lim2) {
            Fmt.print("\nThe $,r term of a[n] is $,d.", count, factors[-1])
            lim2 = lim2 * 10
        }
    }
    i = i + 1
}
Output:
First 100 terms of a[n]:
  1   2   3   2   5   3   7   3   5  11
  3  13   7   5  17   3  19   5   7  11
 23   5  13   7  29   5  31  11  17   7
  3  37  19  13  41   7  43  11   5  23
 47   7   5  17  13  53  11  19  29  59
  5  61  31   7  13  11  67  17  23   7
 71  73  37   5  19  11  13  79  41  83
  7  17  43  29  89   5  13  23  31  47
 19  97   7  11   5 101  17 103   7  53
107 109  11  37 113  19  23  29  13  59

The 1,000th term of a[n] is 109.

The 10,000th term of a[n] is 101.

The 100,000th term of a[n] is 1,693.

The 1,000,000th term of a[n] is 1,202,057.

The 10,000,000th term of a[n] is 1,202,057.
Retrieved from "https://rosettacode.org/wiki/Greatest_prime_dividing_the_n-th_cubefree_number?oldid=361286"