/* Copyright 2012 Kevin Ryde
This file is part of Math-PlanePath.
Math-PlanePath is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the Free
Software Foundation; either version 3, or (at your option) any later
version.
Math-PlanePath is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License along
with Math-PlanePath. If not, see <http://www.gnu.org/licenses/>.
*/
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
typedef unsigned long my_unsigned;
typedef long long my_signed;
#define MY_SIGNED_ABS llabs
#define HYPOT_LIMIT 0x7FFFFFFF
char *
to_base (unsigned long long n, int radix)
{
static char str[256];
static char dstr[256];
int pos = sizeof(str)-1;
do {
int digit = n % radix;
n /= radix;
sprintf (dstr, "[%d]", digit);
int dlen = strlen(dstr);
pos -= dlen;
memcpy (str+pos, dstr, dlen);
} while (n);
return str+pos;
}
int
base_len (unsigned long long n, int radix)
{
int len = 0;
while (n) {
n /= radix;
len++;
}
return len;
}
int
main (void)
{
int realpart, level;
for (realpart = 3; realpart < 10; realpart++) {
int norm = realpart*realpart + 1;
int level_limit = 20;
if (realpart == 2) level_limit = 10;
if (realpart == 3) level_limit = 9;
if (realpart == 4) level_limit = 9;
for (level = 0; level < level_limit; level++) {
unsigned long long min_h = ~0ULL;
my_unsigned min_n = 0;
my_signed min_x = 0;
my_signed min_y = 0;
{
my_unsigned lo = pow(norm, level);
my_unsigned hi = lo * norm;
printf ("%2d lo=%lu hi=%lu\n", level, lo, hi);
my_unsigned n;
for (n = lo; n < hi; n++) {
my_signed x = 0;
my_signed y = 0;
my_signed bx = 1;
my_signed by = 0;
my_unsigned digits = n;
while (digits != 0) {
int digit = digits % norm;
digits /= norm;
x += digit * bx;
y += digit * by;
/* (bx,by) = (bx + i*by)*(i-$realpart) */
my_signed new_bx = bx*-realpart - by;
my_signed new_by = bx + by*-realpart;
bx = new_bx;
by = new_by;
}
unsigned long long abs_x = MY_SIGNED_ABS(x);
unsigned long long abs_y = MY_SIGNED_ABS(y);
if (abs_x > HYPOT_LIMIT
|| abs_y > HYPOT_LIMIT) {
continue;
}
unsigned long long h = abs_x*abs_x + abs_y*abs_y;
/* printf ("%2d %lu %Ld,%Ld %LX\n", level, n, x,y, h); */
if (h < min_h) {
min_h = h;
min_n = n;
min_x = abs_x;
min_y = abs_y;
}
}
}
/* printf ("%lX %Ld,%Ld %s\n", min_n, min_x,min_y, */
/* binary(min_h)); */
printf ("%2d", level);
printf (" %s [%d]", to_base(min_h,norm), base_len(min_h,norm));
printf ("\n");
/* printf ("\n"); */
}
}
return 0;
}