Mysterious logical error - semiprime calculation

[eklavya sharma 2]

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.

//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.

[Smac89]

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)

[eklavya sharma 2]

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]

Here's an idea: how about factoring prod after incrementing count? That way you can easily find the incorrect pair.

[JLBorges]

(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.)

#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).