#include<u.h>
#include<libc.h>
#include<linalg.h>
const float meps = 5.96e-8;
//const double meps = 1.11e-16;
#define MIN(a, b) ((a) < (b) ? (a) : (b))
#define sign(x) (((x) >= 0) - ((x) < 0))
uint MAXITERS = 1000000;
static void *
emalloc(ulong size){
void *v;
v = malloc(size);
if(v == nil)
sysfatal("emalloc: allocating %ld: %r", size);
return v;
}
static void*
erealloc(void *v, ulong sz)
{
v = realloc(v, sz);
if(v == nil)
sysfatal("erealloc: allocating %ld: %r", sz);
return v;
}
static void *
cmalloc(void *p, ulong sz)
{
if(p != nil)
return p;
p = emalloc(sz);
memset(p, 0, sz);
return p;
}
ushort
undilate2(uint x)
{
x = (x | (x >> 1)) & 0x33333333;
x = (x | (x >> 2)) & 0x0F0F0F0F;
x = (x | (x >> 4)) & 0x00FF00FF;
x = (x | (x >> 8)) & 0x0000FFFF;
return x;
}
uint
dilate2(ushort x)
{
uint r = x;
r = (r | (r << 8)) & 0x00FF00FF;
r = (r | (r << 4)) & 0x0F0F0F0F;
r = (r | (r << 2)) & 0x33333333;
r = (r | (r << 1)) & 0x55555555;
return r;
}
/*ushort
undilate3(uint x)
{
x = (x * 0x00015) & 0x0E070381;
x = (x * 0x01041) & 0x0FF80001;
x = (x * 0x40001) & 0x0FFC0000;
return (ushort)(x >> 18);
}
uint
dilate3(ushort x)
{
uint r = x;
r = (r * 0x10001) & 0xFF0000FF;
r = (r * 0x00101) & 0x0F00F00F;
r = (r * 0x00011) & 0xC30C30C3;
r = (r * 0x00005) & 0x49249249;
return r;
}*/
void
fblock(block *p)
{
free(p->data);
free(p);
}
block *
zblock(block *p)
{
memset(p->data, 0, p->sz);
return p;
}
static void
growblock(block *p, ulong sz)
{
if(p->data == nil)
p->allocated = 1;
if(p->sz < sz && p->allocated){
p->data = erealloc(p->data, sz);
p->sz = sz;
}
}
block *
cblock(block *dst, block *src)
{
growblock(dst, src->sz);
memcpy(dst->data, src->data, src->sz);
return dst;
}
vector*
nvector(vector *p, real *data, uint l)
{
ulong sz;
<