Herhangi Matematiksel Şeyler: Euler Projesi 268. Soru

100'den küçük en az dört farklı asal çarpanı bulunan sayıları sayma

100'den küçük en az dört farklı asal çarpanı bulunan 1000'den küçük 23 pozitif tam sayı bulunduğu gösterilebilir.

100'den küçük en az dört farklı asal çarpanı bulunan 10¹⁶'dan küçük kaç pozitif tam sayı vardır?
Cevap: 785478606870985

Mathematica:
let _ =
  let ps = [|2; 3; 5; 7; 11; 13; 17; 19; 23; 29; 31; 37; 41; 43; 47; 53; 59;
      61; 67; 71; 73; 79; 83; 89; 97|] in
  let pk = Array.length ps in
  let n = 10000000000000000 in
  let ss = Array.init (pk + 1) (fun i -> (i - 3) * (i - 2) * (i - 1) / 6) in
  let res = ref 0 in
  let rec bf k l m s =
    if m <= n then (
      if l >= 4
      then res := !res + s * ss.(l) * (n / m);
      for i = k to pk - 1 do
 bf (i + 1) (l + 1) (ps.(i) * m) (-s)
      done;
    )
  in
    bf 0 0 1 1;
    Printf.printf "%d\n" !res

Python:
import itertools, functools, operator, math

primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]

nmax = 10 ** 16
k = 4

def f(n):
    return sum(nmax // functools.reduce(operator.mul, x, 1)
        for x in itertools.combinations(primes, n))

def c(n, m):
    return math.factorial(n) // math.factorial(m) // math.factorial(n - m)

s = 0
for i in range(k, len(primes) + 1):
    ff = f(i)
    if ff == 0:
        break
    s += (-1) ** (i - k) * c(i - 1, k - 1) * ff

print(s)

Pari/gp:
? K=primepi(100) 
time = 0 ms.
%20 = 25
? PK=primes(K) 
time = 0 ms.
%21 = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
? KFC=vector(K);KFC[4]=1;for(i=4+1,K,KFC[i]=1-sum(j=1,i-1,KFC[j]*binomial(i,j)));KFC
time = 4 ms.
%24 = [0, 0, 0, 1, -4, 10, -20, 35, -56, 84, -120, 165, -220, 286, -364, 455, -560, 680, -816, 969, -1140, 1330, -1540, 1771, -2024]
? prime4(bn,k,{m=1},{prd=1},{l=0})=if(prd>bn,0,if(l<4,0,floor(bn/prd)*KFC[l])+sum(i=m,k,prime4(bn,k,i+1,prd*PK[i],l+1)))
time = 0 ms.
? prime4(1000,K) 
time = 40 ms.
%25 = 23
? prime4(10^16,K) 
time = 3mn, 30,906 ms.
%26 = 785478606870985

Haskell:
import ONeillPrimes (primes)

choose :: (Integral a) => a -> a -> a
choose n r = div (product [n-r+1..n]) $ product [1..r]

coeffs :: (Integral a) => [a]
coeffs = zipWith (*) (cycle [-1,1]).map (`choose` 3) $ [0..]

count :: (Integral a) => a -> [a] -> [a] -> a -> a -> a
count n (c:cs) (p:ps) q s
    | p * q > n = s
    | otherwise = let !s1 = count n (c:cs) ps q $ s + c * (div n (p * q))
                  in count n cs ps (p*q) s1
count _ _ _ _ s = s

p268 :: (Integral a) => a -> a -> a
p268 u n = count n cs ps 1 0
    where ps = takeWhile (<u) primes
          cs = take (length ps) coeffs

main :: IO ()
main = print.p268 100 $ 10^16

Go:
from math import factorial
from Primes import Primes
primes = Primes(100)

max = 10**16 - 1

def products(x=1, n=4, i=0):
    
    if x > max:
        pass
    elif n == 1:
        for p in primes[i:]:
            yield x*p
    else:
        for j in range(i, len(primes) - n + 1):
            for p in products(x*primes[j], n - 1, j+1):
                yield p

def f(n):
    
    return factorial(n-1) // (6*factorial(n-4))

def main():
    
    count = 0
    
    for p in products():
        count += max // p
        
    for i in range(5, 15):
        j = 0
        for p in products(n = i):
            j += max // p
        
        if i % 2 == 0:
            count += j * f(i)
        else:
            count -= j * f(i)
        
 print(count)

Java:
public class e268 {
 static int[] primes ={2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97};
 static long limit = 1000000;//((long) 100000000)*((long) 100000000);
 static long prod(String s){
  long ret= 1;
  for(int i=0;i<s.length();i++){
   if((s.charAt(i)-48)==1)
    ret*=(s.charAt(i)-48)*primes[i+(25-s.length())];
   if(ret>= limit) return limit+10;                         
  }
  return ret;
  
 }
 public static void main(String[] args){
  long[] arr = new long[26];
  long time1 = System.currentTimeMillis();
  System.out.println(prod("1111100000000000000000000"));
  for(long i=1;i<33554432;i++){
   String q = Long.toBinaryString(i);
   long r = prod(q);
   if(i%1000000==1)
    System.out.println(i+" "+q+" "+r);
//   if(Long.bitCount(i)==5 && r<limit)
//    System.out.println(i+" "+q+" "+r);
    
   arr[Long.bitCount(i)]+=limit/r;
  }
  
  long tot=0;
  boolean toggle=true;
  for(int i=4;i<=25;i++){
   if(toggle) tot+=arr[i]*((i-3)*(i-2)*(i-1))/6;
   else tot-= arr[i]*((i-3)*(i-2)*(i-1))/6;
   toggle=!toggle;
  }
  System.out.println(tot);
  long time2 = System.currentTimeMillis();
  for(int i=1;i<=25;i++){
   System.out.print(arr[i]+";");
  }
  System.out.println();
  System.out.println((time2-time1)/1000.0);
 }
}

Basic:
Const maxPrime = 100
Const ceiling As Long = 10000000000000000 - 1
Const primesCount As Integer = 4
Dim primes As Integer() = New Integer() {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97}
Dim result As Long = 0

Sub Main()
    Solve(0, 1, -1)
    Console.WriteLine(result)
End Sub

Sub Solve(ByVal level As Integer, ByVal product As Long, ByVal index As Integer)
    If product <= ceiling Then
        Dim delta = level - primesCount
        If delta >= 0 Then result += CLng((1 - (delta Mod 2) * 2) * (delta + 1) * (delta + 2) * (delta + 3) / 6 * Math.Floor(ceiling / product))
        For i = index + 1 To primes.Count - 1
            Solve(level + 1, product * primes(i), i)
        Next
    End If
End Sub

ml:
let primes =
    [    2.;  3.;  5.;  7.; 11.; 13.; 17.; 19.; 23.; 29.;
        31.; 37.; 41.; 43.; 47.; 53.; 59.; 61.; 67.; 71.;
        73.; 79.; 83.; 89.; 97. ]
let factor =
    [| 0.; 0.; 0.; 0.; 1.; -4.; 10.; -20.; 35.; -56.; 84.; -120.; 165.; -220. |]

let main n =
    let rec loop r c acc = function
        | hd :: tl ->
            let g = c *. hd in
            if g > n then acc else
            let acc = acc +. factor.(r) *. floor (n /. g) in
            loop (r+1) g (loop r c acc tl) tl
        | [] -> acc in
    loop 1 1. 0. primes
;;
Printf.printf "Answer: %.f\n" (main (float_of_string Sys.argv.(1)))

Herhangi Matematiksel Şeyler

Sayfalar

7 Temmuz 2019

Euler Projesi 268. Soru