1
0
mirror of https://github.com/zenorogue/hyperrogue.git synced 2024-11-23 21:07:17 +00:00

nilv:: preliminary implementation

This commit is contained in:
Zeno Rogue 2019-08-06 12:00:46 +02:00
parent 3605a7a4e1
commit 4f27b12ca2
19 changed files with 476 additions and 226 deletions

View File

@ -233,7 +233,7 @@ void display_data::set_projection(int ed) {
shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::halfplane3;
if(DIM == 3 && hyperbolic && apply_models && pmodel == mdPerspective)
shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardH3, pers3 = true;
if(DIM == 3 && (euclid || sol) && apply_models && pmodel == mdPerspective)
if(DIM == 3 && (euclid || sol || nil) && apply_models && pmodel == mdPerspective)
shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardR3, pers3 = true;
if(DIM == 3 && apply_models && pmodel == mdSolPerspective)
shaderside_projection = true, glhr::new_shader_projection = glhr::shader_projection::standardSolv, pers3 = true;
@ -334,8 +334,8 @@ void display_data::set_projection(int ed) {
if(pers3) {
glhr::projection_multiply(glhr::frustum(current_display->tanfov, current_display->tanfov * cd->ysize / cd->xsize));
glhr::projection_multiply(glhr::scale(1, -1, -1));
if(solv::local_perspective_used())
glhr::projection_multiply(glhr::tmtogl_transpose(solv::local_perspective));
if(nisot::local_perspective_used())
glhr::projection_multiply(glhr::tmtogl_transpose(nisot::local_perspective));
}
else if(DIM == 3) {
glhr::glmatrix M = glhr::ortho(cd->xsize/current_display->radius/2, -cd->ysize/current_display->radius/2, 1);
@ -344,6 +344,8 @@ void display_data::set_projection(int ed) {
M[2][2] = 2 / (clip_max - clip_min);
M[3][2] = (clip_min + clip_max) / (clip_max - clip_min);
glhr::projection_multiply(M);
if(nisot::local_perspective_used())
glhr::projection_multiply(glhr::tmtogl_transpose(nisot::local_perspective));
}
else {
glhr::projection_multiply(glhr::frustum(cd->xsize / cd->ysize, 1));

View File

@ -1240,7 +1240,7 @@ bool good_for_wall(cell *c) {
}
void buildBigStuff(cell *c, cell *from) {
if(sphere || quotient || sol || (penrose && !binarytiling)) return;
if(sphere || quotient || nonisotropic || (penrose && !binarytiling)) return;
if(chaosmode > 1) return;
bool deepOcean = deep_ocean_at(c, from);

View File

@ -220,6 +220,7 @@ void initcells() {
hrmap* res = callhandlers((hrmap*)nullptr, hooks_newmap);
if(res) currentmap = res;
else if(nonisotropic) currentmap = nisot::new_map();
#if CAP_CRYSTAL
else if(geometry == gCrystal) currentmap = crystal::new_map();
#endif
@ -241,7 +242,6 @@ void initcells() {
else if(sphere) currentmap = new hrmap_spherical;
else if(quotient) currentmap = new quotientspace::hrmap_quotient;
#if CAP_BT
else if(sol) currentmap = solv::new_map();
else if(binarytiling) currentmap = binary::new_map();
#endif
else currentmap = new hrmap_hyperbolic;
@ -407,6 +407,7 @@ int compdist(int dx[]) {
int celldist(cell *c) {
if(fulltorus && WDIM == 2)
return torusmap()->dists[decodeId(c->master)];
if(nil) return DISTANCE_UNKNOWN;
if(euwrap)
return torusconfig::cyldist(decodeId(c->master), 0);
if(masterless)

View File

@ -563,6 +563,7 @@ vector<geometryinfo> ginf = {
{"sol", "none", "Solv geometry", "sol", 8, 3, qBINARY, gcSol, 0x41600, {{7, 5}}, eVariation::pure},
{"kd2", "none", "kite-and-dart", "kd2", 4, 3, qPENROSE, gcEuclid, 0x48000, {{7, 7}}, eVariation::pure},
{"kd3", "none", "kite-and-dart on horospheres", "kd3", 12, 3, qsBP, gcHyperbolic, 0x48200, {{7, 3}}, eVariation::pure},
{"nil", "none", "Nil geometry", "nil", 22, 3, 0, gcNil, 0x48600, {{7, 5}}, eVariation::pure},
//{"product","none", "product space", "product", 7, 3, 0, gcProduct, 0x48400, {{7, 3}}, eVariation::pure},
};

View File

@ -201,10 +201,10 @@ enum eGeometry {
gHoroTris, gHoroRec, gHoroHex,
gField435, gField534,
gBinary4, gSol,
gKiteDart2, gKiteDart3, PROD2(gProduct,)
gKiteDart2, gKiteDart3, gNil, PROD2(gProduct,)
gGUARD};
enum eGeometryClass { gcHyperbolic, gcEuclid, gcSphere, gcSol, PROD(gcProduct) };
enum eGeometryClass { gcHyperbolic, gcEuclid, gcSphere, gcSol, gcNil, PROD(gcProduct) };
enum class eVariation { bitruncated, pure, goldberg, irregular, dual };

View File

@ -465,7 +465,7 @@ void initConfig() {
addsaver(bounded_mine_percentage, "bounded_mine_percentage");
addsaver(solv::geodesic_movement, "solv_geodesic_movement", true);
addsaver(nisot::geodesic_movement, "solv_geodesic_movement", true);
#if CAP_SHMUP
shmup::initConfig();

View File

@ -780,8 +780,8 @@ namespace conformal {
dialog::addBreak(100);
if(sol)
dialog::addBoolItem_action(XLAT("geodesic movement in Sol"), solv::geodesic_movement, 'G');
if(nonisotropic)
dialog::addBoolItem_action(XLAT("geodesic movement in Sol/Nil"), nisot::geodesic_movement, 'G');
dialog::addBoolItem(XLAT("rotation"), do_rotate == 2, 'r');
if(do_rotate == 0) dialog::lastItem().value = XLAT("NEVER");

View File

@ -64,7 +64,7 @@ movedir vectodir(const hyperpoint& P) {
transmatrix U = ggmatrix(cwt.at);
if(GDIM == 3 && WDIM == 2) U = radar_transform * U;
if(solv::local_perspective_used()) U = solv::local_perspective * U;
if(nisot::local_perspective_used()) U = nisot::local_perspective * U;
hyperpoint H = sphereflip * tC0(U);
transmatrix Centered = sphereflip * rgpushxto0(H);

View File

@ -793,6 +793,20 @@ void debug_this() { }
void dqi_poly::draw() {
if(flags & POLY_DEBUG) debug_this();
if(nil && vid.usingGL && pmodel == mdPerspective && (current_display->set_all(global_projection), shaderside_projection)) {
auto npoly = *this;
glcoords.clear();
for(int i=0; i<cnt; i++)
glcoords.push_back(glhr::pointtogl(nisot::inverse_exp(V * glhr::gltopoint( (*tab)[offset+i]), nisot::iTable)));
npoly.offset = 0;
npoly.tab = &glcoords;
npoly.V = Id;
set_width(1);
npoly.gldraw();
return;
}
dynamicval<ld> bs(hr::band_shift, band_shift);
if(!hyperbolic && among(pmodel, mdPolygonal, mdPolynomial)) {
bool any = false;

View File

@ -725,8 +725,8 @@ void showEuclideanMenu() {
}
dialog::addSelItem(XLAT("projection"), current_proj_name(), '1');
dialog::add_action_push(conformal::model_menu);
if(sol)
dialog::addBoolItem_action(XLAT("geodesic movement in Sol"), solv::geodesic_movement, 'G');
if(nonisotropic)
dialog::addBoolItem_action(XLAT("geodesic movement in Sol/Nil"), nisot::geodesic_movement, 'G');
#if CAP_CRYSTAL && MAXMDIM >= 4
crystal::add_crystal_transform('x');
#endif
@ -803,6 +803,10 @@ void showEuclideanMenu() {
dialog::addSelItem(XLAT("Curvature"), XLAT("Sol"), 0);
break;
case gcNil:
dialog::addSelItem(XLAT("Curvature"), XLAT("Nil"), 0);
break;
PROD( case gcProduct:
dialog::addSelItem(XLAT("Curvature"), XLAT("Product"), 0);
break; )

View File

@ -135,7 +135,7 @@ void geometry_information::prepare_basics() {
#endif
#if CAP_BT
if(binarytiling) hexvdist = rhexf = 1, tessf = 1, scalefactor = 1, crossf = hcrossf7;
if(geometry == gHoroRec || penrose || sol) hexvdist = rhexf = .5, tessf = .5, scalefactor = .5, crossf = hcrossf7/2;
if(geometry == gHoroRec || penrose || sol || nil) hexvdist = rhexf = .5, tessf = .5, scalefactor = .5, crossf = hcrossf7/2;
#endif
#if CAP_BT && MAXMDIM >= 4
if(binarytiling) binary::build_tmatrix();

View File

@ -769,13 +769,13 @@ pair<bool, hyperpoint> makeradar(hyperpoint h) {
using namespace hyperpoint_vec;
ld d = hdist0(h);
if(sol && solv::geodesic_movement) {
h = solv::inverse_exp(h, true);
if(sol && nisot::geodesic_movement) {
h = nisot::inverse_exp(h, nisot::iLazy);
ld r = hypot_d(3, h);
if(r < 1) h = h * (atanh(r) / r);
else return {false, h};
}
if(solv::local_perspective_used()) h = solv::local_perspective * h;
if(nisot::local_perspective_used()) h = nisot::local_perspective * h;
if(WDIM == 3) {
if(d >= vid.radarrange) return {false, h};
@ -820,7 +820,7 @@ color_t kind_outline(eItem it) {
transmatrix face_the_player(const transmatrix V) {
if(DIM == 2) return V;
if(sol) return V * cspin(0, 2, ptick(618, 0));
if(nonisotropic) return V * cspin(0, 2, ptick(618, 0));
return rgpushxto0(tC0(V));
}
@ -4429,6 +4429,7 @@ int get_darkval(int d) {
if(binarytiling) return darkval_hbt[d];
if(hyperbolic && S7 == 6) return darkval_e6[d];
if(hyperbolic && S7 == 12) return darkval_s12[d];
if(nil) return ((d % 11) * 3) % 7;
return 0;
}
@ -4504,7 +4505,7 @@ void make_clipping_planes() {
hyperpoint sx = point3(y1 * z2 - y2 * z1, z1 * x2 - z2 * x1, x1 * y2 - x2 * y1);
sx /= hypot_d(3, sx);
sx[3] = 0;
if(solv::local_perspective_used()) sx = inverse(solv::local_perspective) * sx;
if(nisot::local_perspective_used()) sx = inverse(nisot::local_perspective) * sx;
clipping_planes.push_back(sx);
};
ld tx = current_display->tanfov;
@ -5047,7 +5048,7 @@ void drawcell(cell *c, transmatrix V, int spinv, bool mirrored) {
}
if(just_gmatrix) return;
#if MAXMDIM >= 4
if(WDIM == 3 && pmodel == mdPerspective) {
if(WDIM == 3 && pmodel == mdPerspective && !nil) {
using namespace hyperpoint_vec;
hyperpoint H = tC0(V);
for(hyperpoint& cpoint: clipping_planes) if((H|cpoint) < -sin_auto(cgi.corner_bonus)) {
@ -5061,7 +5062,7 @@ void drawcell(cell *c, transmatrix V, int spinv, bool mirrored) {
hyperpoint H = tC0(V);
if(abs(H[0]) <= 2 && abs(H[1]) <= 2 && abs(H[2]) <= 2) ;
else {
hyperpoint H2 = solv::inverse_exp(H, true);
hyperpoint H2 = nisot::inverse_exp(H, nisot::iLazy);
for(hyperpoint& cpoint: clipping_planes) if((H2|cpoint) < -.2) return;
}
noclipped++;
@ -5986,7 +5987,7 @@ void drawcell(cell *c, transmatrix V, int spinv, bool mirrored) {
for(int a=0; a<c->type; a++)
if(c->move(a) && !isWall3(c->move(a), dummy)) {
if(pmodel == mdPerspective && !sphere && !quotient && !penrose && !sol) {
if(pmodel == mdPerspective && !sphere && !quotient && !penrose && !nonisotropic) {
if(a < 4 && among(geometry, gHoroTris, gBinary3) && celldistAlt(c) >= celldistAlt(viewctr.at->c7)) continue;
else if(a < 2 && among(geometry, gHoroRec) && celldistAlt(c) >= celldistAlt(viewctr.at->c7)) continue;
else if(c->move(a)->master->distance > c->master->distance && c->master->distance > viewctr.at->distance && !quotient) continue;
@ -7144,11 +7145,11 @@ void make_actual_view() {
actual_view_transform = solmul(zpush(wall_radar((masterless ? centerover.at : viewctr.at->c7), inverse(View), max)), actual_view_transform * View) * inverse(View);
camera_level = asin_auto(tC0(inverse(actual_view_transform * View))[2]);
}
if(sol) {
if(nonisotropic) {
transmatrix T = actual_view_transform * View;
transmatrix T2 = eupush( tC0(inverse(T)) );
solv::local_perspective = T * T2;
actual_view_transform = inverse(solv::local_perspective) * actual_view_transform;
nisot::local_perspective = T * T2;
actual_view_transform = inverse(nisot::local_perspective) * actual_view_transform;
}
#endif
}
@ -7938,7 +7939,7 @@ cell *viewcenter() {
bool inscreenrange(cell *c) {
if(sphere) return true;
if(euclid) return celldistance(viewcenter(), c) <= get_sightrange_ambush();
if(sol) return gmatrix.count(c);
if(nonisotropic) return gmatrix.count(c);
return heptdistance(viewcenter(), c) <= 8;
}

31
hyper.h
View File

@ -112,8 +112,10 @@ void addMessage(string s, char spamtype = 0);
#define euclid (cgclass == gcEuclid)
#define sphere (cgclass == gcSphere)
#define sol (cgclass == gcSol)
#define nil (cgclass == gcNil)
#define hyperbolic (cgclass == gcHyperbolic)
#define eusol (euclid || sol)
#define nonisotropic (sol || nil)
#define translatable (euclid || nonisotropic)
#define nonorientable (ginf[geometry].flags & qNONORIENTABLE)
#define elliptic (ginf[geometry].flags & qELLIPTIC)
#define quotient (ginf[geometry].flags & qANYQ)
@ -552,7 +554,7 @@ struct gcell {
#define NOBARRIERS 127
#define MODFIXER (2*10090080*17)
#define MAX_EDGE 18
#define MAX_EDGE 22
template<class T> struct walker;
@ -5639,19 +5641,32 @@ namespace kite {
}
#endif
namespace solv {
/* nonisotropic */
namespace nisot {
extern transmatrix local_perspective;
inline bool local_perspective_used() { return sol; }
inline bool local_perspective_used() { return nonisotropic; }
hrmap *new_map();
hyperpoint inverse_exp(hyperpoint h, bool lazy);
transmatrix translate(const hyperpoint h);
bool in_table_range(hyperpoint h);
transmatrix get_solmul(const transmatrix T, const transmatrix V);
transmatrix get_solmul_pt(const transmatrix Position, const transmatrix T);
enum iePrecision { iLazy, iTable };
transmatrix parallel_transport(const transmatrix Position, const transmatrix T);
transmatrix transport_view(const transmatrix T, const transmatrix V);
transmatrix spin_towards(const transmatrix Position, const hyperpoint goal);
hyperpoint inverse_exp(const hyperpoint h, iePrecision p);
}
namespace solv {
extern string solshader;
}
namespace nilv {
extern array<vector<hyperpoint>, 22> facevertices;
void software_renderer(dqi_poly *p);
}
bool in_perspective();
extern int noclipped;

View File

@ -142,11 +142,11 @@ ld edge_of_triangle_with_angles(ld alpha, ld beta, ld gamma) {
// (this is analogous to representing a sphere with points such that x^2+y^2+z^2 == 1)
hyperpoint hpxy(ld x, ld y) {
return hpxyz(x,y, eusol ? 1 : sphere ? sqrt(1-x*x-y*y) : sqrt(1+x*x+y*y));
return hpxyz(x,y, translatable ? 1 : sphere ? sqrt(1-x*x-y*y) : sqrt(1+x*x+y*y));
}
hyperpoint hpxy3(ld x, ld y, ld z) {
return hpxyz3(x,y,z, eusol ? 1 : sphere ? sqrt(1-x*x-y*y-z*z) : sqrt(1+x*x+y*y+z*z));
return hpxyz3(x,y,z, translatable ? 1 : sphere ? sqrt(1-x*x-y*y-z*z) : sqrt(1+x*x+y*y+z*z));
}
// origin of the hyperbolic plane
@ -201,7 +201,7 @@ ld dhypot_d(int d, const hyperpoint& a, const hyperpoint& b) {
}
ld zlevel(const hyperpoint &h) {
if(eusol) return h[GDIM];
if(translatable) return h[GDIM];
else if(sphere) return sqrt(intval(h, Hypc));
else return (h[GDIM] < 0 ? -1 : 1) * sqrt(-intval(h, Hypc));
}
@ -278,20 +278,15 @@ transmatrix eupush(ld x, ld y) {
return T;
}
transmatrix eupush3(ld x, ld y, ld z) {
transmatrix eupush(hyperpoint h) {
if(nonisotropic) return nisot::translate(h);
transmatrix T = Id;
T[0][GDIM] = x;
T[1][GDIM] = y;
if(GDIM == 3) T[2][GDIM] = z;
if(sol) T[0][0] = exp(-z), T[1][1] = exp(+z);
for(int i=0; i<GDIM; i++) T[i][GDIM] = h[i];
return T;
}
transmatrix eupush(hyperpoint h) {
transmatrix T = Id;
for(int i=0; i<GDIM; i++) T[i][GDIM] = h[i];
if(sol) T[0][0] = exp(-h[2]), T[1][1] = exp(+h[2]);
return T;
transmatrix eupush3(ld x, ld y, ld z) {
return eupush(point3(x, y, z));
}
transmatrix euscalezoom(hyperpoint h) {
@ -312,12 +307,8 @@ transmatrix euaffine(hyperpoint h) {
transmatrix cpush(int cid, ld alpha) {
transmatrix T = Id;
if(sol) {
T[cid][GDIM] = alpha;
if(cid == 2)
T[0][0] = exp(-alpha), T[1][1] = exp(+alpha);
return T;
}
if(sol || nil)
return eupush3(cid == 0 ? alpha : 0, cid == 1 ? alpha : 0, cid == 2 ? alpha : 0);
T[GDIM][GDIM] = T[cid][cid] = cos_auto(alpha);
T[cid][GDIM] = sin_auto(alpha);
T[GDIM][cid] = -curvature() * sin_auto(alpha);
@ -339,7 +330,8 @@ bool eqmatrix(transmatrix A, transmatrix B, ld eps) {
// in the 3D space, move the point h orthogonally to the (x,y) plane by z units
hyperpoint orthogonal_move(const hyperpoint& h, ld z) {
if(!hyperbolic) return rgpushxto0(h) * cpush(2, z) * C0;
if(eusol) return hpxy3(h[0], h[1], h[2] + z);
if(nil) return nisot::translate(h) * cpush0(2, z);
if(translatable) return hpxy3(h[0], h[1], h[2] + z);
ld u = 1;
if(h[2]) z += asin_auto(h[2]), u /= acos_auto(z);
u *= cos_auto(z);
@ -473,7 +465,7 @@ transmatrix rpushxto0(const hyperpoint& H) {
}
transmatrix ggpushxto0(const hyperpoint& H, ld co) {
if(eusol) {
if(translatable) {
using namespace hyperpoint_vec;
return eupush(co * H);
}
@ -505,7 +497,7 @@ transmatrix rgpushxto0(const hyperpoint& H) {
// (without using this, imprecision could accumulate)
void fixmatrix(transmatrix& T) {
if(sol) ; // fixed inside solmul
if(nonisotropic) ; // T may be inverse... do not do that
else if(euclid) {
for(int x=0; x<GDIM; x++) for(int y=0; y<=x; y++) {
ld dp = 0;
@ -784,17 +776,19 @@ bool asign(ld y1, ld y2) { return signum(y1) != signum(y2); }
ld xcross(ld x1, ld y1, ld x2, ld y2) { return x1 + (x2 - x1) * y1 / (y1 - y2); }
transmatrix solmul(const transmatrix T, const transmatrix V) {
if(sol) return solv::get_solmul(T, V);
if(nonisotropic) return nisot::transport_view(T, V);
else return T * V;
}
transmatrix solmul_pt(const transmatrix Position, const transmatrix T) {
if(sol) return solv::get_solmul_pt(Position, T);
if(nonisotropic) return nisot::parallel_transport(Position, T);
else return Position * T;
}
transmatrix spin_towards(const transmatrix Position, const hyperpoint goal, int dir, int back) {
transmatrix T = sol ? solv::spin_towards(Position, goal) : rspintox(inverse(Position) * goal);
transmatrix T =
nonisotropic ? nisot::spin_towards(Position, goal) :
rspintox(inverse(Position) * goal);
if(back < 0) T = T * spin(M_PI);
if(dir) T = T * cspin(dir, 0, -M_PI/2);
T = Position * T;

View File

@ -313,7 +313,7 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
switch(pmodel) {
case mdPerspective: {
ld ratio = vid.xres / current_display->tanfov / current_display->radius / 2;
if(solv::local_perspective_used()) H = solv::local_perspective * H;
if(nisot::local_perspective_used()) H = nisot::local_perspective * H;
ret[0] = H[0]/H[2] * ratio;
ret[1] = H[1]/H[2] * ratio;
ret[2] = 1;
@ -321,8 +321,8 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
}
case mdSolPerspective: {
auto S = solv::inverse_exp(H, false);
if(solv::local_perspective_used()) S = solv::local_perspective * S;
auto S = nisot::inverse_exp(H, nisot::iTable);
if(nisot::local_perspective_used()) S = nisot::local_perspective * S;
ld ratio = vid.xres / current_display->tanfov / current_display->radius / 2;
ret[0] = S[0]/S[2] * ratio;
ret[1] = S[1]/S[2] * ratio;
@ -345,7 +345,7 @@ void applymodel(hyperpoint H, hyperpoint& ret) {
}
case mdDisk: {
if(solv::local_perspective_used()) H = solv::local_perspective * H;
if(nisot::local_perspective_used()) H = nisot::local_perspective * H;
ld tz = get_tz(H);
if(!vid.camera_angle) {
ret[0] = H[0] / tz;
@ -886,8 +886,8 @@ transmatrix actualV(const heptspin& hs, const transmatrix& V) {
bool point_behind(hyperpoint h) {
if(sphere) return false;
if(!in_perspective()) return false;
if(pmodel == mdSolPerspective) h = solv::inverse_exp(h, true);
if(solv::local_perspective_used()) h = solv::local_perspective * h;
if(pmodel == mdSolPerspective) h = nisot::inverse_exp(h, nisot::iLazy);
if(nisot::local_perspective_used()) h = nisot::local_perspective * h;
return h[2] < 0;
}
@ -910,6 +910,7 @@ bool invalid_point(const hyperpoint h) {
bool in_smart_range(const transmatrix& T) {
hyperpoint h = tC0(T);
if(invalid_point(h)) return false;
if(nil) return cells_drawn < 2000;
if(pmodel == mdSolPerspective) return solv::in_table_range(h);
hyperpoint h1;
applymodel(h, h1);
@ -1610,7 +1611,7 @@ void queuestraight(hyperpoint X, int style, color_t lc, color_t fc, PPR p) {
void draw_boundary(int w) {
if(w == 1) return;
if(sol) return;
if(nonisotropic || euclid) return;
dynamicval<ld> lw(vid.linewidth, vid.linewidth * vid.multiplier_ring);
@ -1880,7 +1881,11 @@ bool do_draw(cell *c, const transmatrix& T) {
if(WDIM == 3) {
if(cells_drawn > vid.cells_drawn_limit) return false;
if(pmodel == mdSolPerspective) {
if(!solv::in_table_range(tC0(T))) return false;
if(!nisot::in_table_range(tC0(T))) return false;
if(!limited_generation(c)) return false;
}
else if(nil) {
if(!nisot::in_table_range(tC0(T))) return false;
if(!limited_generation(c)) return false;
}
else if(vid.use_smart_range) {

View File

@ -5,20 +5,37 @@
namespace hr {
namespace nisot {
typedef array<float, 3> ptlow;
transmatrix local_perspective;
bool geodesic_movement = true;
transmatrix translate(hyperpoint h) {
transmatrix T = Id;
for(int i=0; i<GDIM; i++) T[i][GDIM] = h[i];
if(sol) {
T[0][0] = exp(-h[2]);
T[1][1] = exp(+h[2]);
}
if(nil)
T[2][1] = h[0];
return T;
}
}
namespace solv {
int PRECX, PRECY, PRECZ;
typedef array<float, 3> ptlow;
vector<ptlow> inverse_exp_table;
vector<nisot::ptlow> inverse_exp_table;
bool table_loaded;
string solfname = "solv-geodesics.dat";
hyperpoint inverse_exp(hyperpoint h);
void load_table() {
if(table_loaded) return;
FILE *f = fopen(solfname.c_str(), "rb");
@ -28,73 +45,18 @@ namespace solv {
fread(&PRECY, 4, 1, f);
fread(&PRECZ, 4, 1, f);
inverse_exp_table.resize(PRECX * PRECY * PRECZ);
fread(&inverse_exp_table[0], sizeof(ptlow) * PRECX * PRECY * PRECZ, 1, f);
fread(&inverse_exp_table[0], sizeof(nisot::ptlow) * PRECX * PRECY * PRECZ, 1, f);
fclose(f);
table_loaded = true;
}
void step2(ld& x, ld &y, ld &z, ld &vx, ld &vy, ld &vz) {
ld ax = -2 * vx * vz;
ld ay = 2 * vy * vz;
ld az = exp(2*z) * vx * vx - exp(-2*z) * vy * vy;
// ld x2 = x + vx/2;
// ld y2 = y + vy/2;
ld z2 = z + vz/2;
ld vx2 = vx + ax/2;
ld vy2 = vy + ay/2;
ld vz2 = vz + az/2;
ld ax2 = -2 * vx2 * vz2;
ld ay2 = 2 * vy2 * vz2;
ld az2 = exp(2*z2) * vx2 * vx2 - exp(-2*z2) * vy2 * vy2;
x += vx + ax/2;
y += vy + ay/2;
z += vz + az/2;
vx += ax2;
vy += ay2;
vz += az2;
}
hyperpoint direct_exp(hyperpoint v, int steps) {
hyperpoint at;
at[0] = 0;
at[1] = 0;
at[2] = 0;
at[3] = 1;
v[0] /= steps;
v[1] /= steps;
v[2] /= steps;
for(int i=0; i<steps; i++) step2(at[0], at[1], at[2], v[0],v[1],v[2]);
return at;
}
void parallel(ld x, ld y, ld z, ld vx, ld vy, ld vz, ld& wx, ld& wy, ld& wz, ld t) {
// ld dax = -dz * gxyz * ay;
ld dwx = -vz * wx - vx * wz;
ld dwy = vz * wy + vy * wz;
ld dwz = vx * wx * exp(2*z) - vy * wy * exp(-2*z);
wx += dwx * t;
wy += dwy * t;
wz += dwz * t;
}
hyperpoint direct_exp(hyperpoint v, int steps, vector<hyperpoint> transported) {
ld x = 0, y = 0, z = 0;
v[0] /= steps;
v[1] /= steps;
v[2] /= steps;
for(int i=0; i<steps; i++) {
for(auto t: transported) parallel(x,y,z, v[0],v[1],v[2], t[0], t[1], t[2], 1);
step2(x,y,z, v[0],v[1],v[2]);
}
return hpxy3(x, y, z);
hyperpoint christoffel(const hyperpoint at, const hyperpoint velocity, const hyperpoint transported) {
return hpxyz3(
-velocity[2] * transported[0] - velocity[0] * transported[2],
velocity[2] * transported[1] + velocity[1] * transported[2],
velocity[0] * transported[0] * exp(2*at[2]) - velocity[1] * transported[1] * exp(-2*at[2]),
0
);
}
ld x_to_ix(ld u) {
@ -111,7 +73,7 @@ namespace solv {
return 0.5 - atan((0.5-x) / y) / M_PI;
}
hyperpoint inverse_exp(hyperpoint h, bool lazy) {
hyperpoint get_inverse_exp(hyperpoint h, bool lazy) {
load_table();
ld ix = h[0] >= 0. ? x_to_ix(h[0]) : x_to_ix(-h[0]);
@ -160,10 +122,6 @@ namespace solv {
return res * atanh(r) / r; */
}
transmatrix local_perspective;
bool geodesic_movement = true;
struct hrmap_sol : hrmap {
hrmap *binary_map;
unordered_map<pair<heptagon*, heptagon*>, heptagon*> at;
@ -299,7 +257,15 @@ namespace solv {
};
hrmap *new_map() { return new hrmap_sol; }
pair<heptagon*,heptagon*> getcoord(heptagon *h) {
return ((hrmap_sol*)currentmap)->coords[h];
}
heptagon *get_at(heptagon *h1, heptagon *h2, bool gen) {
auto m = ((hrmap_sol*)currentmap);
if(!gen && !m->at.count(make_pair(h1, h2))) return nullptr;
return m->get_at(h1, h2);
}
ld solrange_xy = 15;
ld solrange_z = 4;
@ -311,62 +277,6 @@ namespace solv {
return abs(h[0]) < solrange_xy && abs(h[1]) < solrange_xy && abs(h[2]) < solrange_z;
}
void fixmatrix(transmatrix& T) {
transmatrix push = eupush( tC0(T) );
transmatrix push_back = inverse(push);
transmatrix gtl = push_back * T;
{ dynamicval<eGeometry> g(geometry, gSphere); hr::fixmatrix(gtl); }
T = push * gtl;
}
transmatrix spt(transmatrix Pos, transmatrix T) {
solv::fixmatrix(Pos);
hyperpoint h = tC0(T);
h[3] = 0;
h = Pos * h;
int steps = 100;
using namespace hyperpoint_vec;
h /= steps;
for(int i=0; i<steps; i++) {
for(int j=0; j<3; j++)
parallel(Pos[0][3], Pos[1][3], Pos[2][3],
h[0], h[1], h[2],
Pos[0][j], Pos[1][j], Pos[2][j],
+1);
step2(Pos[0][3], Pos[1][3], Pos[2][3], h[0], h[1], h[2]);
}
return Pos;
}
transmatrix get_solmul(const transmatrix T, const transmatrix V) {
if(!geodesic_movement) return V * eupush(inverse(V) * T * V * C0);
return inverse(spt(inverse(V), inverse(T)));
}
transmatrix get_solmul_pt(const transmatrix Position, const transmatrix T) {
if(!geodesic_movement) return Position * T;
return spt(Position, T);
}
transmatrix spin_towards(const transmatrix Position, const hyperpoint goal) {
hyperpoint at = tC0(Position);
transmatrix push_back = inverse(eupush(at));
hyperpoint back_goal = push_back * goal;
back_goal = inverse_exp(back_goal, false);
transmatrix back_Position = push_back * Position;
return rspintox(inverse(back_Position) * back_goal);
}
int approx_distance(heptagon *h1, heptagon *h2) {
auto m = (hrmap_sol*) currentmap;
dynamicval<eGeometry> g(geometry, gBinary4);
@ -412,35 +322,334 @@ namespace solv {
"return res;"
"}";
}
auto sol_config = addHook(hooks_args, 0, [] () {
using namespace arg;
if(argis("-solrange")) {
shift_arg_formula(solrange_xy);
shift_arg_formula(solrange_z);
return 0;
namespace nilv {
hyperpoint christoffel(const hyperpoint Position, const hyperpoint Velocity, const hyperpoint Transported) {
ld x = Position[0];
return point3(
x * Velocity[1] * Transported[1] - 0.5 * (Velocity[1] * Transported[2] + Velocity[2] * Transported[1]),
-.5 * x * (Velocity[1] * Transported[0] + Velocity[0] * Transported[1]) + .5 * (Velocity[2] * Transported[0] + Velocity[0] * Transported[2]),
-.5 * (x*x-1) * (Velocity[1] * Transported[0] + Velocity[0] * Transported[1]) + .5 * x * (Velocity[2] * Transported[0] + Velocity[0] * Transported[2])
);
}
else if(argis("-fsol")) {
shift(); solfname = args();
return 0;
hyperpoint formula_exp(hyperpoint v) {
// copying Modelling Nil-geometry in Euclidean Space with Software Presentation
// v[0] = c cos alpha
// v[1] = c sin alpha
// v[2] = w
if(v[0] == 0 && v[1] == 0) return v;
if(v[2] == 0) return point31(v[0], v[1], v[0] * v[1] / 2);
ld alpha = atan2(v[1], v[0]);
ld w = v[2];
ld c = hypot(v[0], v[1]) / v[2];
return point31(
2 * c * sin(w/2) * cos(w/2 + alpha),
2 * c * sin(w/2) * sin(w/2 + alpha),
w * (1 + (c*c/2) * ((1 - sin(w)/w) + (1-cos(w))/w * sin(w + 2 * alpha)))
);
}
else if(argis("-solglitch")) {
shift_arg_formula(glitch_xy);
shift_arg_formula(glitch_z);
return 0;
hyperpoint get_inverse_exp(hyperpoint h, int iterations) {
ld wmin, wmax;
ld side = h[2] - h[0] * h[1] / 2;
if(hypot_d(2, h) < 1e-6) return h;
else if(side > 1e-6) {
wmin = 0, wmax = 2 * M_PI;
}
else if(side < -1e-6) {
wmin = - 2 * M_PI, wmax = 0;
}
else return point3(h[0], h[1], 0);
ld alpha_total = atan2(h[1], h[0]);
ld cmul;
if(abs(h[0]) > abs(h[1]))
cmul = h[0] / 2 / cos(alpha_total);
else
cmul = h[1] / 2 / sin(alpha_total);
ld sa = sin(2 * alpha_total);
for(int it=0;; it++) {
ld w = (wmin + wmax) / 2;
ld c = cmul / sin(w/2);
ld z = w * (1 + (c*c/2) * ((1 - sin(w)/w) + (1-cos(w))/w * sa));
if(it == iterations) {
ld alpha = alpha_total - w/2;
return point3(c * w * cos(alpha), c * w * sin(alpha), w);
}
if(h[2] > z) wmin = w;
else wmax = w;
}
}
else if(argis("-solgeo")) {
geodesic_movement = true;
pmodel = mdSolPerspective;
return 0;
struct mvec : array<int, 3> {
mvec() { }
mvec(int x, int y, int z) {
auto& a = *this;
a[0] = x; a[1] = y; a[2] = z;
}
mvec inverse() {
auto& a = *this;
return mvec(-a[0], -a[1], -a[2]+a[1] * a[0]);
}
mvec operator * (const mvec b) {
auto& a = *this;
return mvec(a[0] + b[0], a[1] + b[1], a[2] + b[2] + a[0] * b[1]);
}
};
static const mvec mvec_zero = mvec(0, 0, 0);
hyperpoint mvec_to_point(mvec m) { return hpxy3(m[0], m[1], m[2]); }
array<mvec, 22> movevectors = {mvec(0,-1,-1), mvec(1,-1,-1), mvec(-1,0,-1), mvec(0,0,-1), mvec(1,0,-1), mvec(-1,1,-1), mvec(0,1,-1), mvec(-1,-1,0), mvec(0,-1,0), mvec(-1,0,0), mvec(1,1,0), mvec(0,1,1), mvec(-1,1,0), mvec(1,0,1), mvec(0,0,1), mvec(-1,0,1), mvec(1,-1,0), mvec(0,-1,1), mvec(1,1,1), mvec(0,1,0), mvec(1,0,0), mvec(-1,-1,1), };
array<vector<hyperpoint>, 22> facevertices = {{
{point31(0.434375,-0.5,-0.489063), point31(0.414062,-0.417969,-0.578125), point31(0.003125,-0.49375,-0.49375), point31(0.414062,-0.574219,-0.414062), },
{point31(0.434375,-0.5,-0.489063), point31(0.490451,-0.560764,-0.292535), point31(0.414062,-0.574219,-0.414062), },
{point31(-0.49375,0,-0.49375), point31(-0.575,-0.415625,-0.170313), point31(-0.500977,-0.433594,-0.270508), point31(-0.414062,-0.414062,-0.414062), },
{point31(0.414062,-0.417969,-0.578125), point31(0.498047,0.00585938,-0.492188), point31(0.414062,0.414062,-0.414062), point31(-0.0178571,0.5,-0.497768), point31(-0.414062,0.417969,-0.578125), point31(-0.49375,0,-0.49375), point31(-0.414062,-0.414062,-0.414062), point31(0.003125,-0.49375,-0.49375), },
{point31(0.498047,0.00585938,-0.492188), point31(0.575,0.415625,-0.170313), point31(0.503906,0.434896,-0.265625), point31(0.414062,0.414062,-0.414062), },
{point31(-0.414062,0.574219,-0.414062), point31(-0.489258,0.560547,-0.293945), point31(-0.435268,0.504464,-0.485491), },
{point31(-0.0178571,0.5,-0.497768), point31(-0.414062,0.574219,-0.414062), point31(-0.435268,0.504464,-0.485491), point31(-0.414062,0.417969,-0.578125), },
{point31(-0.575,-0.415625,-0.170313), point31(-0.558594,-0.491211,-0.0136719), point31(-0.500977,-0.433594,-0.270508), },
{point31(0.003125,-0.49375,-0.49375), point31(-0.414062,-0.414062,-0.414062), point31(-0.500977,-0.433594,-0.270508), point31(-0.558594,-0.491211,-0.0136719), point31(-0.489258,-0.560547,0.293945), point31(-0.414062,-0.574219,0.414062), point31(-0.00260417,-0.492188,0.497396), point31(0.414062,-0.414062,0.414062), point31(0.503906,-0.434896,0.265625), point31(0.558594,-0.491211,0.0136719), point31(0.490451,-0.560764,-0.292535), point31(0.414062,-0.574219,-0.414062), },
{point31(-0.575,-0.415625,-0.170313), point31(-0.558594,-0.491211,-0.0136719), point31(-0.489258,-0.560547,0.293945), point31(-0.434375,-0.5,0.489063), point31(-0.414062,-0.417969,0.578125), point31(-0.49375,0.0015625,0.495312), point31(-0.575,0.415625,0.170313), point31(-0.558594,0.492969,0.0109375), point31(-0.489258,0.560547,-0.293945), point31(-0.435268,0.504464,-0.485491), point31(-0.414062,0.417969,-0.578125), point31(-0.49375,0,-0.49375), },
{point31(0.575,0.415625,-0.170313), point31(0.558594,0.492969,-0.0109375), point31(0.503906,0.434896,-0.265625), },
{point31(0.435268,0.504464,0.485491), point31(0.414062,0.574219,0.414062), point31(0.0146484,0.499023,0.499023), point31(0.414062,0.417969,0.578125), },
{point31(-0.558594,0.492969,0.0109375), point31(-0.500977,0.433594,0.270508), point31(-0.575,0.415625,0.170313), },
{point31(0.503906,-0.434896,0.265625), point31(0.575,-0.415625,0.170313), point31(0.497396,-0.00390625,0.494792), point31(0.414062,-0.414062,0.414062), },
{point31(-0.00260417,-0.492188,0.497396), point31(-0.414062,-0.417969,0.578125), point31(-0.49375,0.0015625,0.495312), point31(-0.414062,0.414062,0.414062), point31(0.0146484,0.499023,0.499023), point31(0.414062,0.417969,0.578125), point31(0.497396,-0.00390625,0.494792), point31(0.414062,-0.414062,0.414062), },
{point31(-0.49375,0.0015625,0.495312), point31(-0.414062,0.414062,0.414062), point31(-0.500977,0.433594,0.270508), point31(-0.575,0.415625,0.170313), },
{point31(0.558594,-0.491211,0.0136719), point31(0.575,-0.415625,0.170313), point31(0.503906,-0.434896,0.265625), },
{point31(-0.414062,-0.574219,0.414062), point31(-0.434375,-0.5,0.489063), point31(-0.414062,-0.417969,0.578125), point31(-0.00260417,-0.492188,0.497396), },
{point31(0.490451,0.560764,0.292535), point31(0.414062,0.574219,0.414062), point31(0.435268,0.504464,0.485491), },
{point31(-0.0178571,0.5,-0.497768), point31(-0.414062,0.574219,-0.414062), point31(-0.489258,0.560547,-0.293945), point31(-0.558594,0.492969,0.0109375), point31(-0.500977,0.433594,0.270508), point31(-0.414062,0.414062,0.414062), point31(0.0146484,0.499023,0.499023), point31(0.414062,0.574219,0.414062), point31(0.490451,0.560764,0.292535), point31(0.558594,0.492969,-0.0109375), point31(0.503906,0.434896,-0.265625), point31(0.414062,0.414062,-0.414062), },
{point31(0.414062,-0.417969,-0.578125), point31(0.498047,0.00585938,-0.492188), point31(0.575,0.415625,-0.170313), point31(0.558594,0.492969,-0.0109375), point31(0.490451,0.560764,0.292535), point31(0.435268,0.504464,0.485491), point31(0.414062,0.417969,0.578125), point31(0.497396,-0.00390625,0.494792), point31(0.575,-0.415625,0.170313), point31(0.558594,-0.491211,0.0136719), point31(0.490451,-0.560764,-0.292535), point31(0.434375,-0.5,-0.489063), },
{point31(-0.414062,-0.574219,0.414062), point31(-0.434375,-0.5,0.489063), point31(-0.489258,-0.560547,0.293945), },
}};
struct hrmap_nil : hrmap {
unordered_map<mvec, heptagon*> at;
unordered_map<heptagon*, mvec> coords;
heptagon *getOrigin() override { return get_at(mvec_zero); }
heptagon *get_at(mvec c) {
auto& h = at[c];
if(h) return h;
h = tailored_alloc<heptagon> (S7);
h->c7 = newCell(S7, h);
coords[h] = c;
h->dm4 = 0;
h->zebraval = c[0];
h->emeraldval = c[1];
h->fieldval = c[2];
h->cdata = NULL;
h->alt = NULL;
return h;
}
heptagon *create_step(heptagon *parent, int d) override {
auto p = coords[parent];
auto q = p * movevectors[d];
auto child = get_at(q);
parent->c.connect(d, child, (d + 11) % 22, false);
return child;
}
transmatrix adjmatrix(int i) {
return nisot::translate(mvec_to_point(movevectors[i]));
}
virtual transmatrix relative_matrix(heptagon *h2, heptagon *h1) override {
return nisot::translate(mvec_to_point(coords[h1].inverse() * coords[h2]));
}
void draw() override {
dq::visited.clear();
dq::enqueue(viewctr.at, cview());
while(!dq::drawqueue.empty()) {
auto& p = dq::drawqueue.front();
heptagon *h = get<0>(p);
transmatrix V = get<1>(p);
dq::drawqueue.pop();
cell *c = h->c7;
if(!do_draw(c, V)) continue;
drawcell(c, V, 0, false);
if(0) for(int t=0; t<c->type; t++) {
if(!c->move(t)) continue;
dynamicval<color_t> g(poly_outline, darkena((0x142968*t) & 0xFFFFFF, 0, 255) );
queuepoly(V, cgi.shWireframe3D[t], 0);
}
for(int i=0; i<S7; i++) {
// note: need do cmove before c.spin
heptagon *h1 = h->cmove(i);
dq::enqueue(h1, V * adjmatrix(i));
}
}
}
};
}
namespace nisot {
hyperpoint christoffel(const hyperpoint at, const hyperpoint velocity, const hyperpoint transported) {
if(sol) return solv::christoffel(at, velocity, transported);
else if(nil) return nilv::christoffel(at, velocity, transported);
else return point3(0, 0, 0);
}
else if(argis("-solnogeo")) {
geodesic_movement = false;
pmodel = mdPerspective;
return 0;
bool in_table_range(hyperpoint h) {
if(sol) return solv::in_table_range(h);
return true;
}
return 1;
});
hyperpoint inverse_exp(const hyperpoint h, iePrecision p) {
if(sol) return solv::get_inverse_exp(h, p == iLazy);
if(nil) return nilv::get_inverse_exp(h, 10);
return point3(h[0], h[1], h[2]);
}
void geodesic_step(hyperpoint& at, hyperpoint& velocity) {
using namespace hyperpoint_vec;
auto acc = christoffel(at, velocity, velocity);
auto at2 = at + velocity / 2;
auto velocity2 = velocity + acc / 2;
auto acc2 = christoffel(at2, velocity2, velocity2);
at = at + velocity + acc2 / 2;
velocity = velocity + acc;
}
hyperpoint direct_exp(hyperpoint v, int steps) {
using namespace hyperpoint_vec;
hyperpoint at = point31(0, 0, 0);
v /= steps;
v[3] = 0;
for(int i=0; i<steps; i++) geodesic_step(at, v);
return at;
}
transmatrix transpose(transmatrix T) {
transmatrix result;
for(int i=0; i<MDIM; i++)
for(int j=0; j<MDIM; j++)
result[j][i] = T[i][j];
return result;
}
transmatrix parallel_transport_bare(transmatrix Pos, transmatrix T) {
hyperpoint h = tC0(T);
h[3] = 0;
h = Pos * h;
int steps = 100;
using namespace hyperpoint_vec;
h /= steps;
auto tPos = transpose(Pos);
for(int i=0; i<steps; i++) {
for(int j=0; j<3; j++)
tPos[j] += christoffel(tPos[3], h, tPos[j]);
geodesic_step(tPos[3], h);
}
return transpose(tPos);
}
void fixmatrix(transmatrix& T) {
transmatrix push = eupush( tC0(T) );
transmatrix push_back = inverse(push);
transmatrix gtl = push_back * T;
{ dynamicval<eGeometry> g(geometry, gSphere); hr::fixmatrix(gtl); }
T = push * gtl;
}
transmatrix parallel_transport(const transmatrix Position, const transmatrix T) {
if(!geodesic_movement) return Position * T;
auto P = Position;
nisot::fixmatrix(P);
return parallel_transport_bare(P, T);
}
transmatrix transport_view(const transmatrix T, const transmatrix V) {
if(!geodesic_movement) return V * eupush(inverse(V) * T * V * C0);
return inverse(parallel_transport(inverse(V), inverse(T)));
}
transmatrix spin_towards(const transmatrix Position, const hyperpoint goal) {
hyperpoint at = tC0(Position);
transmatrix push_back = inverse(translate(at));
hyperpoint back_goal = push_back * goal;
back_goal = inverse_exp(back_goal, iTable);
transmatrix back_Position = push_back * Position;
return rspintox(inverse(back_Position) * back_goal);
}
hrmap *new_map() {
if(sol) return new solv::hrmap_sol;
if(nil) return new nilv::hrmap_nil;
return NULL;
}
auto config = addHook(hooks_args, 0, [] () {
using namespace arg;
if(argis("-solrange")) {
shift_arg_formula(solv::solrange_xy);
shift_arg_formula(solv::solrange_z);
return 0;
}
else if(argis("-fsol")) {
shift(); solv::solfname = args();
return 0;
}
else if(argis("-solglitch")) {
shift_arg_formula(solv::glitch_xy);
shift_arg_formula(solv::glitch_z);
return 0;
}
else if(argis("-solgeo")) {
geodesic_movement = true;
pmodel = mdSolPerspective;
return 0;
}
else if(argis("-solnogeo")) {
geodesic_movement = false;
pmodel = mdPerspective;
return 0;
}
return 1;
});
}

View File

@ -367,7 +367,7 @@ int fieldval_uniq(cell *c) {
auto p = cell_to_pair(c);
return gmod(p.first * torusconfig::dx + p.second * torusconfig::dy, torusconfig::qty);
}
else if(binarytiling || archimedean) return 0;
else if(binarytiling || archimedean || nil) return 0;
else if(&currfp == &fp_invalid) return 0;
else if(WDIM == 3) return c->master->fieldval;
else if(ctof(c) || NONSTDVAR) return c->master->fieldval/S7;

View File

@ -909,7 +909,7 @@ void geometry_information::create_wall3d() {
}
}
if(DIM == 3 && !euclid && !binarytiling) {
if(DIM == 3 && !euclid && !binarytiling && !nil) {
reg3::generate();
int facesize = isize(reg3::cellshape) / S7;
for(int w=0; w<S7; w++) {
@ -935,6 +935,10 @@ void geometry_information::create_wall3d() {
make_wall(7, {pt(-1,00,+1), pt(+1,00,+1), pt(+1,-1,+1), pt(-1,-1,+1)});
}
if(geometry == gNil) {
for(int i=0; i<S7; i++) make_wall(i, nilv::facevertices[i]);
}
if(penrose) {
auto kv = kite::make_walls();
for(auto& v: kv.first) for(auto& h: v) {

View File

@ -948,7 +948,7 @@ bincode acd_bin(ld x) {
bincode get_bincode(hyperpoint h) {
switch(ginf[gwhere].cclass) {
case gcEuclid: case gcSol: PROD( case gcProduct: )
case gcEuclid: case gcSol: case gcNil: PROD( case gcProduct: )
return acd_bin(h[0]) + acd_bin(h[1]) * sY + acd_bin(h[2]) * sZ;
case gcHyperbolic:
return acd_bin(hypot_d(3, h));