Author Topic: Testing if a rational is in the cantor set (Read 6250 times)

ruler501 · « **on:** February 03, 2014, 05:05:58 pm »

Is there a fast way to test for membership in the cantor set. My only idea is to test out to some threshold(25 digits or so) and if there are no ones(in ternary) before that to assume it is in the cantor set. This might work, but it is likely to include more rationals than it should(though q will always be under 3^20). Is there a better(preferably fast) way to do this? My code currently spends 97% of it's time in this function(though it is called ~2 million times) and any improvements will give me massive performance boosts.

This is my current code(THRESH is 70 which should be more precision than unsigned long long provides)

Code: [Select]

bool testCantor(unsigned long long p, unsigned long long q)
{
    unsigned long double rem = ((long double)p)/q, div = ((long double)1)/3;
    int test;
    unsigned long long mult = 3;
    for(int i=0; i<THRESH; i++)
    {
        test = rem * mult;
        if(test == 1)
        {
             if(rem > div) return false;
             else return true;
        }
        rem -= test*div;
        mult *= 3;
        div /= 3;
    }
    return true;
}

calcdude84se · « **Reply #1 on:** February 04, 2014, 02:07:40 pm »

I'm not sure why you'd want to use floating-point. You can just use integer math.
The key to the following code is that, if you write out p/q in ternary, because it is rational it will eventually repeat.
This could probably use some optimization, but here's some code that should work (except for the noted caveats):

Code: [Select]

bool testCantor(unsigned long long p, unsigned long long q) // 0 <= p <= q
{
  //First check all the non-repeating digits
  while(q % 3 == 0)
  {
    q /= 3;
    if(p / q == 1)
      return p == q;
    else
      p %= q;
  }
  unsigned long long p_start = p;
  //Now test all the repeating digits once
  do
  {
    unsigned long long p3 = p * 3; // Note: if q's current value is over one-third of the maximum value of an unsigned long long, overflow errors will be a problem. You'll need to fix this somehow. This is the line that could overflow.
    if(p3 / q == 1)
      return false;
    else
      p = p3 % q;
  } while(p != p_start);
  return true;
}

I'm not sure how fast it is, but it should return the correct answer.
Edit: Changed the first return false; to return p == q;, because numbers of the form 0.xxxxx1 (or 0.xxxxx02222...) are in the Cantor set.

ruler501 · « **Reply #2 on:** February 04, 2014, 02:35:58 pm »

It runs slightly faster than mine did and is much more accurate. One problem though your check for p/q==1 needs to return p==q not false because if it terminates with a 1 it is in the Cantor set (change the 1 to a 0 and add repeating 2 after it)

calcdude84se · « **Reply #3 on:** February 04, 2014, 03:19:50 pm »

Ah, right, I forgot about the edges being slightly tricky. Was that the only problem? I'll fix my code for posterity.

ruler501 · « **Reply #4 on:** February 04, 2014, 03:32:29 pm »

That's the only problem I found. I'll post my whole program later
Edit: also my code will never give it more than 3^20 and unsigned long long shouldn't overflow till 3^40

EDIT2: Here's the whole program. It is counting all rationals with denominators(in simplest form) between 3^n and 3^(n+1) in the cantor set

Code: [Select]

#include <iostream>
#include <algorithm>
#include <thread>
#include <cmath>

#define THREADS 8

using namespace std;

bool testCantor(unsigned long long p, unsigned long long q){
   while(q % 3 == 0){
      q /= 3;
      if(p/q == 1) return p == q;
      p %= q;
   }
   unsigned long long p_start = p;
   do{
      unsigned long long p3 = p * 3;
      if(p3/q == 1) return false;
      p = p3 % q;
   } while(p != p_start);
   return true;
}

void genRational(unsigned long long minNum, unsigned long long maxNum, int jump, unsigned long long *result){
    for(unsigned long long q = minNum; q <= maxNum; q += jump)
        for(unsigned long long p = 1; p < q/3; p++)
            if(p % 3 != 0 && testCantor(p, q) && __gcd(p,q) == 1)
                for(unsigned long long i = p; i < q/3; i *= 3)
                        (*result) += 2 * (__gcd(i,q) == 1);
}

int main(int argc, char* argv[]){
    int n = argc;
    unsigned long long minNum = pow(3, n), *tempLongLong = nullptr, maxNum = pow(3, n+1), result = 0;
    vector<thread*> usedThreads;
    vector<unsigned long long*> results;
    thread* tempThread = nullptr;
    for(int i = 0; i < THREADS; i++){
        tempLongLong = new unsigned long long(0);
        results.push_back(tempLongLong);
        tempThread = new thread(genRational, minNum + i, maxNum, THREADS, tempLongLong);
        usedThreads.push_back(tempThread);
    }
    for(auto i = usedThreads.begin(); i != usedThreads.end(); i++) (*i)->join();
    for(auto i = results.begin(); i != results.end(); i++) result += **i;
    cout << n << " " << result << endl;
}

EDIT: Used an additional property for halving the search space(if a number x is in the cantor set then 1-x is also in the set). I haven't tested the code yet(using windows which stubbornly won't compile if I'm using threads)

EDIT2: Have it down to only checking 2/9 of the rationals(though depending on rounding may miss 1, but that should barely affect results)
also as a note to calcdude84se your else statements in the function didn't actually do anything(the if returned so execution stopped there if it was true), so I removed them. Not sure if that actually boosts performance though.

calcdude84se · « **Reply #5 on:** February 05, 2014, 11:22:42 am »

Quote from: ruler501 on February 04, 2014, 03:32:29 pm

EDIT2: Have it down to only checking 2/9 of the rationals(though depending on rounding may miss 1, but that should barely affect results)

To speak of possible speed-ups, if you have the same denominator (ignoring powers of 3) a lot, you might be able to get a performance boost by caching the number of times the loop that checks the repeating digits goes around and using that number to count the loop instead. (Not sure if it's significant or worthwhile, though; what it would enable is some loop unrolling)

Quote

also as a note to calcdude84se your else statements in the function didn't actually do anything(the if returned so execution stopped there if it was true), so I removed them. Not sure if that actually boosts performance though.

It's more of a stylistic thing whether to use an else with if()return (though I'm not terribly consistent about my usage. I try to pick what's clearer, but it's not always apparent. I might have gotten it backwards this time.), and I would hope that the compiler generates the same code in each case.

ruler501 · « **Reply #6 on:** February 05, 2014, 01:40:01 pm »

I do have the same denominator a lot(order of 2*3^12 to 2 * 3^13 times in the current run), but I can't cache all of them because that would be 2*3^15 different values(stored in 8 bytes each) which would be 219mb of data(actually this might be reasonable) I might try to find a way to just pass it in so it doesn't have to be stored.
I'll do that later today and speed test it to see if there is any significant benefit.

ruler501 · « **Reply #7 on:** February 12, 2014, 12:28:15 am »

I've now managed to get an OpenCL version of it running, but it seems to run a lot slower than it does on the processor. Is this just not a good problem for GPU's or is there a bigger problem.

Code: [Select]

#define __CL_ENABLE_EXCEPTIONS
#include <CL/cl.hpp>
#include <functional>
#include <ctime>
#include <iostream>
#include <fstream>
#include <exception>
#include <cstdlib>
#include <vector>
#include <thread>
#include <cmath>
#include <string>
#include <algorithm>
#include <thread>
#include <cmath>
#include <sstream>

#define SUCCESS 0
#define FAILURE 1
#define EXPECTED_FAILURE 2

const int NUM_ELEMENTS = 512;

int convertToString(const char *filename, std::string& s)
{
    size_t size;
    char*  str;

    // create a file stream object by filename
    std::fstream f(filename, (std::fstream::in | std::fstream::binary));


    if(!f.is_open())
    {
     	return FAILURE;
    }
    else
    {
        size_t fileSize;
        f.seekg(0, std::fstream::end);
        size = fileSize = (size_t)f.tellg();
        f.seekg(0, std::fstream::beg);

        str = new char[size+1];
        if(!str)
        {
            f.close();
            return FAILURE;
        }

        f.read(str, fileSize);
        f.close();
        str[size] = '\0';

        s = str;
        delete[] str;
        return SUCCESS;
    }
}

void printOutput(unsigned long long start, unsigned long long *values){
    for(unsigned int i = 0; i < NUM_ELEMENTS; i++)
       if (values[i] != 0)
            std::cout << start+i << ',' << values[i] << std::endl;
}

void newList(unsigned long long start, unsigned long long *dataList){
    for(int i=0; i < NUM_ELEMENTS; ++i)
        dataList[i] = start + i;
}

using namespace cl;

Kernel kernelA;
Context context;
CommandQueue queue;

int init() {
    cl_int status = 0;
	const char* buildOption ="-x clc++ ";
	std::vector<Platform> platforms;
	status = Platform::get(&platforms);
	if (status != CL_SUCCESS){
		std::cout<<"Error: Getting platforms!"<<std::endl;
		return FAILURE;
	}
	std::vector<cl::Platform>::iterator iter;
	for(iter = platforms.begin(); iter != platforms.end(); ++iter)
		if(!strcmp((*iter).getInfo<CL_PLATFORM_VENDOR>().c_str(), "Advanced Micro Devices, Inc."))
            break;
	cl_context_properties cps[3] = {CL_CONTEXT_PLATFORM, (cl_context_properties)(*iter) (), 0};
	bool gpuNotFound = false;
	try{
		context = cl::Context(CL_DEVICE_TYPE_GPU, cps, NULL, NULL, &status);
	}
	catch(std::exception e){
		gpuNotFound = true;
	}
	if(gpuNotFound){
		std::cout<<"GPU not found, falling back to CPU!"<<std::endl;
		context = cl::Context(CL_DEVICE_TYPE_CPU, cps, NULL, NULL, &status);
		if (status != CL_SUCCESS){
			std::cout<<"Error: Creating context!"<<std::endl;
			return FAILURE;
		}
	}
	try{
		std::vector<Device> devices = context.getInfo<CL_CONTEXT_DEVICES>();
		queue = CommandQueue(context, devices[0]);
		std::ifstream sourceFile("Rationals.cl");
		std::string sourceCode(
			std::istreambuf_iterator<char>(sourceFile),
			(std::istreambuf_iterator<char>()));
		Program::Sources source(1, std::make_pair(sourceCode.c_str(), sourceCode.length()+1));
		Program program = Program(context, source);
		program.build(devices, buildOption);
		kernelA = Kernel(program, "countRationals");
	}catch(cl::Error e){
		std::cout << "Line "<< __LINE__<<": Error in "<<e.what() <<std::endl;
		return FAILURE;
	}
	return SUCCESS;
}

int execute(unsigned long long* inputList, unsigned long long* outputList) {
    try{
		Buffer inputBuffer = Buffer(context, CL_MEM_READ_WRITE, NUM_ELEMENTS * sizeof(unsigned long long));
		Buffer outputBuffer = Buffer(context, CL_MEM_READ_WRITE, NUM_ELEMENTS * sizeof(unsigned long long));
		queue.enqueueWriteBuffer(inputBuffer, CL_TRUE, 0, NUM_ELEMENTS * sizeof(unsigned long long), inputList);
		kernelA.setArg(0, inputBuffer);
		kernelA.setArg(1, outputBuffer);
		NDRange global(NUM_ELEMENTS);
		NDRange local(NUM_ELEMENTS/2);
		queue.enqueueNDRangeKernel(kernelA, NullRange, global, local);
		queue.enqueueReadBuffer(outputBuffer, CL_TRUE, 0, NUM_ELEMENTS * sizeof(unsigned long long), outputList);
	}catch(cl::Error e){
		std::cout << "Line "<< __LINE__<<": Error in "<<e.what() <<std::endl;
		return FAILURE;
	}
    return SUCCESS;
}

using namespace std;

int main(int argc, char* argv[]){
    unsigned long long minNum, maxNum;
    if (argc == 2){
        minNum = pow(3, atoi(argv[1]));
        maxNum = pow(3, atoi(argv[1]) + 1);
    }
    else if (argc == 3){
        minNum = pow(3, atoi(argv[1]));
        maxNum = pow(3, atoi(argv[2]));
    }
    else if (argc == 4){
        minNum = pow(3, atoi(argv[1]));
        maxNum = pow(3, atoi(argv[2]));
    }
    else return -1;
    unsigned long long *q = nullptr, *result = nullptr, *old = nullptr;
    thread workThread, outThread;
    q = new unsigned long long[NUM_ELEMENTS];
    newList(minNum, q);
    result = new unsigned long long[NUM_ELEMENTS];
    init();
    workThread = thread(execute, q, result);
    workThread.join();
    for(unsigned long long i = minNum + NUM_ELEMENTS; i < maxNum; i += NUM_ELEMENTS){
        old = result;
        result = new unsigned long long[NUM_ELEMENTS];
        newList(i, q);
        workThread = thread(execute, q, result);
        outThread = thread(printOutput, i, old);
        workThread.join();
        outThread.join();
        delete[] old;
        old = nullptr;
    }
    delete[] q;
    delete[] result;
    return 0;
}

With this kernel code

Code: [Select]

bool testCantor(unsigned long p, unsigned long q){
	while(q % 3 == 0){
		q /= 3;
		if (p/q == 1) return p==q;
		p %= q;
	}
	unsigned long p_start = p;
	do{
		unsigned long p3 = p * 3;
		if(p3/q == 1) return false;
		p = p3 % q;
	} while(p != p_start);
	return true;
}

bool coprime(unsigned long a, unsigned long b){
    while (a != b){
        if (a > b) a = a - b;
        else b = b - a;
    }
    return a == 1;
}

__kernel
void countRationals(__global unsigned long *input, __global unsigned long *output){
    int gid = get_global_id(0);
    unsigned long q = input[gid], p = 1;
    output[gid] = 0;
    for(p = 1; p <= q/3; p++){
        if(p % 3 != 0 && testCantor(p, q))
            for(unsigned long i = p; i <= q/3; i *= 3)
                if(coprime(i,q))
                    output[gid] += 2;
    }
}

Author Topic: Testing if a rational is in the cantor set (Read 6250 times)

ruler501

Testing if a rational is in the cantor set

calcdude84se

Re: Testing if a rational is in the cantor set

ruler501

Re: Testing if a rational is in the cantor set

calcdude84se

Re: Testing if a rational is in the cantor set

ruler501

Re: Testing if a rational is in the cantor set

calcdude84se

Re: Testing if a rational is in the cantor set

ruler501

Re: Testing if a rational is in the cantor set

ruler501

Re: Testing if a rational is in the cantor set