Mysterious logical error - semiprime calculation

eklavya sharma 2 (170)
I've been trying to solve Project Euler problem 123 - semiprimes

A composite is a number containing at least 2 prime factors. For example, 15 = 3 * 5; 9 = 3 * 3; 12 = 2 * 2 * 3.

There are 10 composites below 30 containing precisely 2, not necessarily distinct, prime factors: 4, 6, 9, 10, 14, 15, 21, 22, 25, 26.

How many composite integers, n < 10^8, have precisely 2, not necessarily distinct, prime factors?


At first I tried to make a very efficient algorithm, but it gave wrong answer. So I made a brute force algorithm to generate answer for small values of n and then use them to cross check my efficient algorithm and find bugs. I was pretty sure my brute force algorithm would work for all n, but sadly it didn't.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
 
//Semiprimes
#include<iostream>
#include<fstream>
#include<algorithm>
#include<eku\math\sqrt_uint.h>

using std::cout;
using std::cin;
using std::endl;

void init_pbf(unsigned* pbf,unsigned size)
//fills pbf with primes from a file containing 10^8 primes
{
	std::ifstream ifile("E:\\PND\\1e+8.pb32",std::ios::binary);
	if(!ifile.is_open()){std::cout<<"Err";return;}
	ifile.read((char*)pbf,size*sizeof(unsigned));
}

int main()
{
	unsigned i,j,n,size,sqrt_n;
	unsigned long long count,prod;
	unsigned* pbf;
//	unsigned* lb;
	
	cout<<"Enter pbf size: ";
	cin>>size;
//Enter 10^8
	pbf=new unsigned[size];
	init_pbf(pbf,size);
	while(true)
	{
		cout<<"Enter n: ";
		cin>>n;
		sqrt_n=eku::sqrt_uint<unsigned>(n).first;
		count=0;
		/*for(i=0;pbf[i]<=sqrt_n;++i)
		{
			lb=std::upper_bound(pbf,pbf+size,n/pbf[i])-1;
			count+=lb-pbf-i+1;
		}*/
		for(i=0;pbf[i]<=sqrt_n;++i)for(j=i;pbf[j]<=n;++j)
		{
			prod=pbf[i];
			prod*=pbf[j];
			if(prod>=n)break;
			++count;
		}
		cout<<count<<endl;
	}
	return 0;
}

Then I found this OEIS page: http://oeis.org/A066265

This sadly revealed the answer, but I used the values to cross check my brute force approach. It works for n = 10, 100, 1000, 10^4, 10^5, 10^6, 10^7, but not for 10^8. The answer my program displays is one more than the correct answer. I have no idea why this is happening.

The file <eku\math\sqrt_uint.h> is my own and contains a function to calculate square root of an integer. The answer is rounded down if the parameter is not a perfect square.

The file containing the list of primes has never betrayed me in any Project Euler problem before, so there is very little chance that the file is wrong.
Last edited on
Smac89 (1516)
This might be a brute-force approach, but I don`t see any other efficient method of doing this except using hashmaps:

All you need to solve this is all the prime numbers below 10,000 (9973 and below). Store all the values 4-10^8 in a bool array and mark each one as true.

Now using 2 for-loops, go through all the primes you have and as you mark the product of the 2 for-loops as false, keep counting. When your for-loop is done, you have your answer if you print your counter variable.

E: The double for-loops should not scare you as there are only 1229 primes below 10000 so this code should run quite faster than O(1229^2)
Last edited on
eklavya sharma 2 (170)
Well I have devised an O(NlogN) algorithm, but that's not my focus right now. What I'm curious about is that my brute force algorithm doesn't seem to work for n =10^8 and larger n.
helios (11788)
Here's an idea: how about factoring prod after incrementing count? That way you can easily find the incorrect pair.
JLBorges (5443)
(10^8 is a small number. By the prime number theorem, the number of primes less than or equal to 10^8 / 2 is about 2.8 million; we can easily hold them in memory.)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
 
#include <type_traits>
#include <vector>
#include <iostream>

constexpr auto MAX = 100000000 ; // 10^8
using int_type = std::remove_const< decltype(MAX) >::type ;

std::vector<int_type> prime_numbers()
{
    const int_type SZ = MAX / 2 + 1 ;
    std::vector<bool> sieve( SZ, true ) ;
    for( int_type i = 2 ; i < SZ ; ++i ) if( sieve[i] )
       for( int_type j = i+i ; j < SZ ; j += i ) sieve[j] = false ;

    std::vector<int_type> result ;
    for( int_type i = 2 ; i < SZ ; ++i ) if( sieve[i] ) result.push_back(i) ;
    return result ;
}

int main()
{
    auto primes = prime_numbers() ;

    int_type count = 0 ;

    for( auto p : primes )
    {
        for( auto p2 : primes )
        {
            if( p*p2 >= MAX ) break ;
            if( p2 >= p ) ++count ;
        }
    }

    std::cout << count << '\n' ; // 17427258
}

http://ideone.com/I4ibdE

EDIT: Note: The logic is unchanged from the code in the original post; but generates the prime numbers instead of reading them from a file. (And deduces an appropriate unsigned integral type capable of holding numbers up to 10^8).
Last edited on
Topic archived. No new replies allowed.

C++

Forum