(05-07-2021 10:35 AM)Albert Chan Wrote: Since test for bucket n/2 is same as matching gcd = n/2, we can remove it.
Hopefully, this is my final version, by removing stuff from recurse5() / puzzle5().
Nah, there's more to remove.
Odd/even is covered by whether the GCD is a multiple of 2 or not. So the odd/even test can go.
As hinted at in my answer to your private message, as well as in my previous code, the mod4 test on non-mod4 bases can be worked into the GCD partitioning as well. And as I found out after writing the code, it can also be generalized to higher powers of 2, more specifically to the smallest power of 2 that doesn't divide the base anymore.
I'm doubling the base for the purpose of GCD calls - that covers this mod4 case and similar cases with higher powers of 2, but it needs an adjustment after creating the buckets: swap the mod4 one (or whatever power of 2) with the mod2 one (or next-lower power of 2 in the general case), that accounts for the offset of half the mod between qualifying digits and the mod. Exception: if the base is a power of 2 itself, all GCDs are powers of 2 as well, and the additional bucket is outside of our range of digits, so no swapping is needed. The next-lower power of 2 would also be the base itself, which isn't a valid bucket either because the digit range goes to just 1 below that, so depending on how your languages handles out-of-bounds accesses, you could just swap anyway.
To illustrate:
- On base 10:
- - creating buckets with GCD(20, i): 1->{1, 3, 7, 9}, 2->{2, 6}, 4->{4, 8}, 5->{5}
- - highest power of 2 dividing the base is 2, so this extra bucket is mod4 -> the buckets to swap are 2 and 4
- - result: digits allowed in positions {1, 3, 7, 9} are {1, 3, 7, 9}, digits allowed in positions {2, 6} are {4, 8}, digits allowed in positions {4, 8} are {2, 6}, and position 5 has to be 5, of course.
- On base 12:
- - creating buckets with GCD(24, i): 1->{1, 5, 7, 11}, 2->{2, 6, 10}, 4->{4}, 6->{6}, 8->{8}
- - highest power of 2 dividing the base is 4, so this extra bucket is mod8 -> the buckets to swap are 4 and 8
- - result: digits allowed in positions {1, 5, 7, 11} are {1, 5, 7, 11}, digits allowed in positions {2, 6, 10} are {2, 6, 10}, position 4 has to be 8, position 6 has to be 6, and position 8 has to be 4.
- On base 8:
- - creating buckets with GCD(16, i): 1->{1, 3, 5, 7}, 2->{2, 6}, 4->{4} (note: these are the same buckets as with GCD(8, i), doubling the base has no effect on powers of 2)
- - highest power of 2 dividing the base is 8, so this extra bucket would be mod16 -> don't swap to avoid index-out-of-bounds
- - result: digits allowed in positions {1, 3, 5, 7} are {1, 3, 5, 7}, digits allowed in positions {2, 6} are {2, 6}, and position 4 has to be 4.
---
The mentioned private message was an inquiry regarding what turned out to be bugs in my programs. The first is a simple uninitialized variable in nearinflen_uint_iszero, where the declaration of i needs to be turned into a definition to 0. (Why didn't GCC warn me about use of an uninitialized variable?) The function is only used in printing anyway, so the calculations wouldn't be affected by the bug. Either way, it's fixed in the code listing below.
The more serious bug is an incorrectly patched-in mod4 detection, which resulted in no solutions found on bases 6, 10, and 14. (Silly me for not testing the patch sufficiently!) This bug was actually present not only in the C version but also in the patched SysRPL one, since they both use the same algorithm. The explanation above is for the corrected version; they previously only doubled up the base, and did so only if it wasn't divisible by 4. The generalization removes this divisibility-by-4 check, and the bugfix swaps the buckets. Corrected code for both languages follows.
Code:
#include <assert.h>
#include <malloc.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
typedef unsigned int singledata_t;
typedef unsigned long long doubledata_t;
static_assert (sizeof(singledata_t) * 2 == sizeof(doubledata_t));
typedef struct {
unsigned int len;
singledata_t data[];
} nearinflen_uint_t;
#define NEARINFLEN_UINT_MALLOC(len) ((nearinflen_uint_t *) \
malloc(sizeof(nearinflen_uint_t) + sizeof(singledata_t) * (len)))
#define NEARINFLEN_UINT_SHIFT (sizeof(singledata_t) * 8)
static inline nearinflen_uint_t *nearinflen_uint_zero() {
nearinflen_uint_t *n = NEARINFLEN_UINT_MALLOC(1);
assert(n != NULL);
n->len = 1;
n->data[0] = 0;
return n;
}
static inline bool nearinflen_uint_iszero(nearinflen_uint_t *n) {
int i = 0;
for (; i < n->len; ++i)
if (n->data[i] != 0) return false;
return true;
}
static inline void nearinflen_uint_release(nearinflen_uint_t *n) {
free(n);
}
static nearinflen_uint_t *nearinflen_uint_add(nearinflen_uint_t *n, singledata_t m) {
nearinflen_uint_t *res = NEARINFLEN_UINT_MALLOC(n->len + 1);
assert(res != NULL);
unsigned int i = 0;
doubledata_t carry = m;
for (; i < n->len && carry != 0; ++i) {
carry += n->data[i];
res->data[i] = (singledata_t) carry;
carry >>= NEARINFLEN_UINT_SHIFT;
}
if (carry != 0) {
res->data[n->len] = (singledata_t) carry;
res->len = n->len + 1;
} else {
res->len = n->len;
memcpy(res->data + i, n->data + i, sizeof(singledata_t) * (n->len - i));
}
return res;
}
static nearinflen_uint_t *nearinflen_uint_mul(nearinflen_uint_t *n, singledata_t m) {
nearinflen_uint_t *res = NEARINFLEN_UINT_MALLOC(n->len + 1);
assert(res != NULL);
unsigned int i = 0;
doubledata_t carry = 0;
for (; i < n->len; ++i) {
carry += ((doubledata_t) n->data[i]) * m;
res->data[i] = (singledata_t) carry;
carry >>= NEARINFLEN_UINT_SHIFT;
}
if (carry != 0) {
res->data[n->len] = (singledata_t) carry;
res->len = n->len + 1;
} else res->len = n->len;
return res;
}
static singledata_t nearinflen_uint_div2(nearinflen_uint_t *n, singledata_t m,
nearinflen_uint_t **res) {
if (res != NULL) {
*res = NEARINFLEN_UINT_MALLOC(n->len);
assert(*res != NULL);
memset((*res)->data, 0, sizeof(singledata_t) * n->len);
}
unsigned int i = n->len - 1;
doubledata_t rem = n->data[i];
if (m > rem) {
if (i == 0) {
if (res != NULL) (*res)->len = 1;
return n->data[0];
}
--i;
rem = (rem << NEARINFLEN_UINT_SHIFT) | n->data[i];
}
if (res != NULL) (*res)->len = i + 1;
doubledata_t m_shifted = m, bit = 1;
for (; rem >= m_shifted; bit <<= 1, m_shifted <<= 1);
while (bit != 1 || i != 0) {
if (bit == 1) {
bit <<= NEARINFLEN_UINT_SHIFT;
m_shifted <<= NEARINFLEN_UINT_SHIFT;
--i;
rem = (rem << NEARINFLEN_UINT_SHIFT) | n->data[i];
}
bit >>= 1;
m_shifted >>= 1;
if (rem >= m_shifted) {
rem -= m_shifted;
if (res != NULL) (*res)->data[i] |= bit;
}
}
return (singledata_t) rem;
}
typedef unsigned long long bitset_t;
typedef struct output_digit_t output_digit_t;
struct output_digit_t {
singledata_t digit;
output_digit_t *next;
};
singledata_t base;
bitset_t *valid_digits;
static void try_digits(nearinflen_uint_t *n, singledata_t i, bitset_t not_taken) {
bitset_t valid = not_taken & valid_digits[i];
singledata_t digit = i - nearinflen_uint_div2(n, i, NULL);
for (; digit < base; digit += i) {
bitset_t digit_bitset = 1 << (digit - 1);
if ((valid & digit_bitset) == 0) continue;
nearinflen_uint_t *tmp = nearinflen_uint_add(n, digit);
if (i + 1 != base) {
nearinflen_uint_t *tmp2 = nearinflen_uint_mul(tmp, base);
nearinflen_uint_release(tmp);
try_digits(tmp2, i + 1, not_taken & ~digit_bitset);
continue;
}
/* this destroys digit and n, but we cannot have
two results that only differ in the last digit anyway,
so print and break out */
nearinflen_uint_release(n);
digit = nearinflen_uint_div2(tmp, base, &n);
bool last;
output_digit_t *out = NULL, *prev;
do {
last = nearinflen_uint_iszero(n);
prev = out;
out = malloc(sizeof(output_digit_t));
assert(out != NULL);
out->next = prev;
out->digit = digit;
nearinflen_uint_release(tmp);
tmp = n;
digit = nearinflen_uint_div2(tmp, base, &n);
} while (!last);
nearinflen_uint_release(tmp);
do {
last = out->next == NULL;
printf("%u%c", out->digit, last ? '\n' : ';');
prev = out;
out = out->next;
free(prev);
} while (!last);
break;
}
nearinflen_uint_release(n);
}
static inline singledata_t gcd(singledata_t big, singledata_t small) {
do {
singledata_t tmp = big % small;
big = small;
small = tmp;
} while (small != 0);
return big;
}
int main(int argc, char **argv) {
if (argc < 2) {
printf("Usage: %s <num>\n", argv[0]);
return 1;
}
base = (singledata_t) atoi(argv[1]);
if (base < 2 || base > 65) {
printf("Invalid input, must be an integer between 2 and 65\n");
return 1;
}
if ((base & 1) != 0) return 0;
valid_digits = calloc(base, sizeof(bitset_t));
assert(valid_digits != NULL);
singledata_t *gcd_by_pos = calloc(base, sizeof(singledata_t));
assert(gcd_by_pos != NULL);
singledata_t i = 1;
bitset_t bit = 1;
for (; i < base; ++i, bit <<= 1) {
singledata_t tmp = gcd(i, base * 2);
gcd_by_pos[i] = tmp;
valid_digits[tmp] |= bit;
}
singledata_t highest_2power_gcd = 1;
while (base % highest_2power_gcd == 0) highest_2power_gcd *= 2;
if (highest_2power_gcd < base) {
bitset_t tmp = valid_digits[highest_2power_gcd];
valid_digits[highest_2power_gcd] = valid_digits[highest_2power_gcd / 2];
valid_digits[highest_2power_gcd / 2] = tmp;
}
for (i = 1; i < base; ++i) valid_digits[i] = valid_digits[gcd_by_pos[i]];
free(gcd_by_pos);
nearinflen_uint_t *zero = nearinflen_uint_zero();
try_digits(zero, 1, ~((bitset_t) 0));
free(valid_digits);
return 0;
}
Code:
::
CK1NoBlame FPTR2 ^CK1Z
DUP FPTR2 ^Z># DUP #2-
BINT63 #> caseSIZEERR
DUPONE #AND #0<> case2DROP
WORDSIZE 1LAMBIND
ERRSET ::
BINT64 dostws BINT1 #>HXS
OVER ONE_DO
DUP ONE{}N 4UNROLL
DUP4UNROLL bitSL
LOOP
DROP' ::
OVER #2* #3+ GETLAM
4PICK bitAND UNROT
3GETLAM FPTR2 ^QMul 1GETLAM
3PICK 3PICKOVER FPTR2 ^#>Z
FPTR2 ^Mod FPTR2 ^Z># #-
DO
INDEX@ #2* #2+ GETLAM
4PICKOVER bitAND
OVER HXS==HXS %0=
ITE_DROP ::
bitNOT 5PICK bitAND
3PICK#1+_ 3PICK INDEX@
FPTR2 ^#>Z FPTR2 ^QAdd
OVER 1GETLAM #<>case
2GETEVAL
ROTDROPSWAP
#2* #2* #3- UNROLL
;
OVER +LOOP
4DROP
;
SWAP' NULLLAM OVER #2* #1+
NDUPN DOBIND
1GETLAM ONE_DO
1GETLAM #2* INDEX@
BEGIN
SWAPOVER #/ DROP
#0=UNTIL
DROP #2* #3+ DUP GETLAM
INNERCOMP INDEX@ ROTOVER
#2* #2+ GETLAM bitOR
ROT#1+ {}N SWAP PUTLAM
LOOP
BINT1
BEGIN
#2*
1GETLAM OVER #/ DROP
#0<> UNTIL
1GETLAM OVER#<
ITE_DROP ::
DUP #2* #3+
SWAPOVER GETLAM INNERCOMP
get1 #3+
DUP GETLAM INNERDUP
#4+PICK SWAPROT
OVER #4+ UNPICK_
{}N SWAP PUTLAM
{}N SWAP PUTLAM
;
1GETLAM ONE_DO
INDEX@ #2* #3+ GETLAM
DUPTYPEHSTR?
ITE_DROP ::
INNERCOMP ONE_DO
DUPROT #2* #3+ PUTLAM
LOOP
DROP
;
LOOP
BINT0 #>HXS bitNOT BINT1 Z0_
2GETEVAL
ABND
; ERRTRAP ::
1GETLAM dostws ERRJMP
;
1GETABND dostws
;
Of course the bug invalidates the timing and results I got from that iteration of the algorithm. 24 now takes 111.7_s on the 50g, 22 takes 379.3_s. On the laptop, 60 is still rescued by the large number of small buckets and processes in about a second. Out of the numbers below 60, only 58 takes a significant amount of time (still running after an hour away). I don't have accurate measurements, but the other ones are done in less than five minutes, often significantly less.
Edit: another hour and a half later, I noticed the 58 run must have finished while I wasn't looking. With that, we're looking at no solutions on bases from 16 up to 60.
Search
Member List
Calendar
Help
