Ah, direct quadratic interpolation!

Sometimes I just have to slap myself for missing something like that.

Printable View

- 04-30-2012iMalc
Ah, direct quadratic interpolation!

Sometimes I just have to slap myself for missing something like that. - 04-30-2012oogabooga
I notice a bug in my code in post#15. Since it's using unsigned int (uint), it could overflow and not be noticed with a simple i < n comparison (leading to an infinite loop). I'm not sure of the best way to fix that, but perhaps:

Code:`int isPentagonal(uint n) {`

uint i, delta, last_i = 0;

for (i = 1, delta = 4;

i < n && i > last_i;

last_i = i, i += delta, delta += 3)

;

return i == n;

}

- 04-30-2012oogabooga
I decided to implement this and realized that it's pretty simple. The following seems to work fairly well, up to 2 to the 64 (about 18 x 10 to the 18). The calculation of p is divided into even and odd cases of n to avoid overflow.

Code:`#include <stdio.h>`

#include <stdlib.h>

#include <limits.h>

#include <stdbool.h>

typedef unsigned long long ULL;

// Number of pentagonal numbers up to 2 to the 64

#define PENT_N_MAX 3506826112

bool isPentagonal(ULL x) {

ULL p, n, lo = 1, hi = PENT_N_MAX;

while (lo < hi) {

n = (lo + hi) / 2;

if (n % 2 == 0) p = n / 2 * (3 * n - 1);

else p = (3 * n - 1) / 2 * n;

if (x < p) hi = n;

else if (x > p) lo = n + 1;

else return true;

}

return false;

}

int main(void) {

char line[128];

ULL x;

for (;;) {

printf("Enter number: ");

fgets(line, sizeof line, stdin);

x = ULLONG_MAX;

sscanf(line, "%llu", &x);

if (x == 0) break;

if (x != ULLONG_MAX)

printf("%s\n", isPentagonal(x)?"yes":"no");

}

return 0;

}

- 04-30-2012phantomotap
O_o

I like that, but I have to say that your style is revolting.

Your style is revolting. ^_^;

Seriously though, you might get a smudge more performance by removing the multiple of two related branch. You have a binary search so assuming that rounding loses one you still only have at most one extra branch per number instead of per loop.

Soma - 04-30-2012oogabooga
Ouch! I would appreciate a rewrite both for style and the removal of the even/odd branch.

- 04-30-2012phantomotapQuote:

for style

I don't like the use of spacing. It is repellant to my eyes. I especially don't like the use of vertical space doing the job of creating logical blocks of code because that is exactly the job done by braces. I read through it twice thinking it was wrong because of that style, but I didn't figure it really was because I've seen a lot of your posts.

If, for whatever reason, you want to see examples of my code, you'll find dozens of samples if you search for my name.

Quote:

the removal branch

I'm literally just talking about removing the branch and having the same mathematics always performed.

[Edit]

Sorry. I got cut off.

You would certainly need to do division last.

Code:`p = (n * (3 * n - 1)) / 2;`

[/Edit]

Soma - 04-30-2012oogabooga
I see. Fair enough. The fact that it made you do a double-take is evidence enough for me against it. I don't usually space it like that anyway.

As for removing the branch, leaving just one or the other formula without the branch doesn't seem to work.

The basic formula for the n'th pentagonal number is

p(n) = (3*n*n - n) / 2

which works just fine, unless it overflows, which it does for high n. The only solution there being to lower the highest allowable n.

I split it into either (3*n-1) * (n/2) or (3*n-1)/2 * n to avoid overflow. Note that if the n/2 was off by one, the calculated pentagonal number would be off by (3*n-1). So, since either n is even or (3*n-1) is even, I divide the even factor by 2 using the branch.

At any rate, no biggie. And I hope the following code is less repellant. :)

Code:`// pent.c: Determine if the entered number is pentagonal.`

// Formula for the nth pentagonal number is p(n) = (3*n*n - n) / 2

#include <stdio.h>

#include <stdlib.h>

#include <limits.h>

#include <stdbool.h>

typedef unsigned long long ULL;

// Number of pentagonal numbers up to 2 to the 64

#define N_MAX 3506826112

bool isPentagonal(ULL x) {

ULL lo = 1, hi = N_MAX;

while (lo < hi) {

// index of pentagonal number

ULL n = (lo + hi) / 2;

ULL p; // pentagonal number p(n)

if (n % 2 == 0)

p = n/2 * (3*n - 1);

else

p = (3*n - 1)/2 * n;

if (x < p)

hi = n;

else if (x > p)

lo = n + 1;

else

return true;

}

return false;

}

int main(void) {

char line[128];

ULL x;

for (;;) {

printf("Enter number: ");

fgets(line, sizeof line, stdin);

x = ULLONG_MAX;

sscanf(line, "%llu", &x);

if (x == 0)

break;

if (x != ULLONG_MAX)

printf("%s\n", isPentagonal(x)?"yes":"no");

}

return 0;

}

- 04-30-2012phantomotapQuote:

At any rate, no biggie. And I hope the following code is less repellant.

Quote:

The only solution there being to lower the highest allowable n.

*shrug*

Anyway, as I said in one edit, you'd need to reduce the margin of error to just one instead of `(3*n-1)' requiring an extra iteration and an extra layer in the process for very high numbers. Well, you could as you say assume a smaller domain. (I also noted this in a lost edit.)

You have, cleverly, coded a binary search that performs the same operations for every value in the domain. I'm, ultimately, talking about moving some of the branching of the inner most loop out allowing you to assume a lot more than you now can for the lower half of the domain. The higher value half of the domain would be handled by the same code you have now. It would increase the size of the function a little but reduce the number of "tight" branches significantly. This would traditionally help with AMD and hurt a little with Intel depending on the side of the code in the utility set using those inner most loops.

*shrug*

Well, I did say "you might", but I should have noted that it clear that is was just a hunch meaning that it just as likely would fail.

Soma - 05-01-2012phantomotap
Well, actually, you can pretty much always eliminate branches related to mathematics by introducing more mathematics.

After optimization it is sometimes even faster than a simple branching version depending on the processor.

Here it will not be, but because you said my hunch was wrong I'll post it anyway.*****

*pbblbtllbtlt*

Code:`ULL nd2 = n/2;`

ULL nm2 = n%2;

ULL n3p1 = (3*n - 1);

p1 = (nd2 * n3p1) + ((nm2 * n3p1) >> 1);

*****I didn't actually test this either. ^_^;