circle_fitting

Untitled7.ipynb

{
  "cells": [
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "from clifford.g3c import *\nfrom clifford.tools.g3c import *\nfrom clifford.tools.g3c.object_fitting import *\nfrom pyganja import *\nimport numba",
      "execution_count": 1,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "@numba.njit\ndef val_fit_circle(point_list):\n    \"\"\"\n    Performs Leo Dorsts circle fitting technique\n    \"\"\"\n    # Check if there are just 3 points\n    if point_list.shape[0] == 3:\n        best_obj = point_list[0, :]\n        for i in range(1, 3):\n            best_obj = omt_func(best_obj, point_list[i, :])\n        return val_normalised(best_obj)\n    # Loop over our points and construct the matrix\n    accumulator_matrix = np.zeros((32, 32))\n    for i in range(point_list.shape[0]):\n        # Get the point as a left gmt matrix\n        P_i_l = get_left_gmt_matrix(point_list[i, :])\n        # Multiply and add\n        accumulator_matrix += P_i_l @ mask0 @ P_i_l\n    accumulator_matrix = accumulator_matrix @ mask1\n    # Find the eigenvalues of the matrix\n    e_vals, e_vecs = np.linalg.eig(accumulator_matrix)\n#     print(e_vals)\n    # Find the smallest and second smallest non negative eigenvalues\n    min_eval = np.inf\n    min_eval2 = np.inf\n    min_eval_index = -1\n    min_eval_index2 = -1\n    for i in range(len(e_vals)):\n        this_e_val = e_vals[i]\n        if this_e_val > 0:\n            if this_e_val < min_eval:\n                min_eval2 = min_eval\n                min_eval = this_e_val\n                min_eval_index2 = min_eval_index\n                min_eval_index = i\n            else:\n                if this_e_val < min_eval2:\n                    min_eval2 = this_e_val\n                    min_eval_index2 = i\n#     print(min_eval_index, min_eval)\n#     print(min_eval_index2, min_eval2)\n    best_sphere = val_normalised(mask1@np.real(e_vecs[:, min_eval_index]))\n    second_best_sphere = val_normalised(mask1@np.real(e_vecs[:, min_eval_index2]))\n    best_circle = val_normalised(mask3@dual_func(omt_func(best_sphere, second_best_sphere)))\n    return best_circle\n\n\ndef fit_circle(point_list):\n    \"\"\"\n    Performs Leo Dorsts circle fitting technique\n    \"\"\"\n    return layout.MultiVector(value=val_fit_circle(np.array([p.value for p in point_list])))\n",
      "execution_count": 2,
      "outputs": []
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "# Test the algorithms"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "noise = 0.1\nnpnts = 10\n\ntrue_circle = random_circle()\npoint_list = project_points_to_circle([random_conformal_point() for i in range(npnts)], true_circle)\npoint_list = [up(down(P) + noise * random_euc_mv()) for P in point_list]\nsphere = fit_sphere(point_list)\nplane = fit_plane(point_list)\ncircle = meet(sphere,plane).normal()\nprint(true_circle.normal())\nprint(circle.normal())\n\ngs = GanjaScene()\ngs.add_objects(point_list, color=Color.BLACK, static=True)\ngs.add_objects([true_circle], color=Color.GREEN)\ngs.add_objects([circle], color=Color.RED)\ngs.add_objects([fit_circle(point_list)], color=Color.BLUE)\ndraw(gs, scale=0.05)",
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "text": "(0.58049^e123) + (12.14244^e124) + (12.18716^e125) + (0.21685^e134) + (0.22322^e135) + (0.11658^e145) - (6.69146^e234) - (6.70045^e235) + (0.32736^e245) + (0.07009^e345)\n(0.6578^e123) + (12.16264^e124) + (12.20989^e125) - (0.02278^e134) - (0.01596^e135) + (0.12766^e145) - (7.29606^e234) - (7.30503^e235) + (0.35821^e245) + (0.07591^e345)\n",
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": "<IPython.core.display.Javascript object>",
            "application/javascript": "/** Ganja.js - Geometric Algebra - Not Just Algebra. \n  * @author Enki\n  * @link   https://github.com/enkimute/ganja.js\n  */\n\n/*********************************************************************************************************************/\n// \n// Ganja.js is an Algebra generator for javascript. It generates a wide variety of Algebra's and supports operator\n// overloading, algebraic literals and a variety of graphing options.\n//\n// Ganja.js is designed with prototyping and educational purposes in mind. Clean mathematical syntax is the primary\n// target.\n//\n// Ganja.js exports only one function called *Algebra*. This function is used to generate Algebra classes. (say complex\n// numbers, minkowski or 3D CGA). The returned class can be used to create, add, multiply etc, but also to upgrade\n// javascript functions with algebraic literals, operator overloading, vectors, matrices and much more.\n//\n// As a simple example, multiplying two complex numbers 3+2i and 1+4i could be done like this :\n//\n//  var complex = Algebra(0,1);\n//  var a = new complex([3,2]);\n//  var b = new complex([1,3]); \n//  var result = a.Mul(b);\n// \n// But the same can be written using operator overloading and algebraic literals. (where scientific notation with\n// lowercase e is overloaded to directly specify generators (e1, e2, e12, ...))\n//\n//   var result = Algebra(0,1,()=>(3+2e1)*(1+4e1));\n//\n// Please see github for user documentation and examples. \n//\n/*********************************************************************************************************************/\n\n// Documentation below is for implementors. I'll assume you know about Clifford Algebra's, grades, its products, etc ..\n// I'll also assume you are familiar with ES6. My style may feel a bith mathematical, advise is to read slow. \n\n(function (name, context, definition) {\n  if (typeof module != 'undefined' && module.exports) module.exports = definition();\n  else if (typeof define == 'function' && define.amd) define(name, definition);\n  else context[name] = definition();\n}('Algebra', this, function () {\n\n/** The Algebra class generator. Possible calling signatures : \n  *   Algebra([func])                      => algebra with no dimensions, i.e. R. Optional function for the translator.\n  *   Algebra(p,[func])                    => 'p' positive dimensions and an optional function to pass to the translator.\n  *   Algebra(p,q,[func])                  => 'p' positive and 'q' negative dimensions and optional function.\n  *   Algebra(p,q,r,[func])                => 'p' positive, 'q' negative and 'r' zero dimensions and optional function.\n  *   Algebra({                            => for custom basis, cayley, mixing, etc pass in an object as first parameter.\n  *     [p:p],                             => optional 'p' for # of positive dimensions\n  *     [q:q],                             => optional 'q' for # of negative dimensions\n  *     [r:r],                             => optional 'r' for # of zero dimensions\n  *     [metric:array],                    => alternative for p,q,r. e.g. ([1,1,1,-1] for spacetime)\n  *     [basis:array],                     => array of strings with basis names. (e.g. ['1','e1','e2','e12'])\n  *     [Cayley:Cayley],                   => optional custom Cayley table (strings). (e.g. [['1','e1'],['e1','-1']])                            \n  *     [mix:boolean],                     => Allows mixing of various algebras. (for space efficiency).\n  *     [graded:boolean],                  => Use a graded algebra implementation. (automatic for +6D)\n  *     [baseType:Float32Array]            => optional basetype to use. (only for flat generator)\n  *   },[func])                            => optional function for the translator.\n **/  \n  return function Algebra(p,q,r) {\n  // Resolve possible calling signatures so we know the numbers for p,q,r. Last argument can always be a function.\n    var fu=arguments[arguments.length-1],options=p; if (options instanceof Object) {\n      q = (p.q || (p.metric && p.metric.filter(x=>x==-1).length))| 0;\n      r = (p.r || (p.metric && p.metric.filter(x=>x==0).length)) | 0;\n      p = p.p === undefined ? (p.metric && p.metric.filter(x=>x==1).length) || 0 : p.p || 0;\n    } else { options={}; p=p|0; r=r|0; q=q|0; };\n\n  // Support for multi-dual-algebras\n    if (p==0 && q==0 && r<0) { r=-r; // Create a dual number algebra if r>1 .. consider more explicit syntax\n      options.basis  = [...Array(r+1)].map((a,i)=>i?'e0'+i:'1');  options.metric = [1,...Array(r)]; options.tot=r+1;\n      options.Cayley = [...Array(r+1)].map((a,i)=>[...Array(r+1)].map((y,j)=>i*j==0?((i+j)?'e0'+(i+j):'1'):'0'));\n    }\n  \n\n  // Calculate the total number of dimensions.\n    var tot = options.tot = (options.tot||(p||0)+(q||0)+(r||0)||(options.basis&&options.basis.length))|0;\n \n  // Unless specified, generate a full set of Clifford basis names. We generate them as an array of strings by starting\n  // from numbers in binary representation and changing the set bits into their relative position.  \n  // Basis names are ordered first per grade, then lexically (not cyclic!). \n  // For 10 or more dimensions all names will be double digits ! 1e01 instead of 1e1 .. \n    var basis=options.basis||[...Array(2**tot)]                                                                                 // => [undefined, undefined, undefined, undefined, undefined, undefined, undefined, undefined]\n              .map((x,xi)=>(((1<<30)+xi).toString(2)).slice(-tot||-1)                                                           // => [\"000\", \"001\", \"010\", \"011\", \"100\", \"101\", \"110\", \"111\"]  (index of array in base 2)\n              .replace(/./g,(a,ai)=>a=='0'?'':String.fromCharCode(66+ai-(r!=0))))                                               // => [\"\", \"3\", \"2\", \"23\", \"1\", \"13\", \"12\", \"123\"] (1 bits replaced with their positions, 0's removed)\n              .sort((a,b)=>(a.toString().length==b.toString().length)?(a>b?1:b>a?-1:0):a.toString().length-b.toString().length) // => [\"\", \"1\", \"2\", \"3\", \"12\", \"13\", \"23\", \"123\"] (sorted numerically)\n              .map(x=>x&&'e'+(x.replace(/./g,x=>('0'+(x.charCodeAt(0)-65)).slice(tot>9?-2:-1) ))||'1')                          // => [\"1\", \"e1\", \"e2\", \"e3\", \"e12\", \"e13\", \"e23\", \"e123\"] (converted to commonly used basis names)\n              \n  // See if the basis names start from 0 or 1, store grade per component and lowest component per grade.             \n    var low=basis.length==1?1:basis[1].match(/\\d+/g)[0]*1,\n        grades=options.grades||basis.map(x=>tot>9?(x.length-1)/2:x.length-1),\n        grade_start=grades.map((a,b,c)=>c[b-1]!=a?b:-1).filter(x=>x+1).concat([basis.length]);\n\n  // String-simplify a concatenation of two basis blades. (and supports custom basis names e.g. e21 instead of e12)      \n  // This is the function that implements e1e1 = +1/-1/0 and e1e2=-e2e1. The brm function creates the remap dictionary.\n    var simplify = (s,p,q,r)=>{\n          var sign=1,c,l,t=[],f=true,ss=s.match(tot>9?/(\\d\\d)/g:/(\\d)/g);if (!ss) return s; s=ss; l=s.length;\n          while (f) { f=false;\n          // implement Ex*Ex = metric.\n            for (var i=0; i<l;) if (s[i]===s[i+1]) { if ((s[i]-low)>=(p+r)) sign*=-1; else if ((s[i]-low)<r) sign=0; i+=2; f=true; } else t.push(s[i++]);\n          // implement Ex*Ey = -Ey*Ex while sorting basis vectors.  \n            for (var i=0; i<t.length-1; i++) if (t[i]>t[i+1]) { c=t[i];t[i]=t[i+1];t[i+1]=c;sign*=-1;f=true; break;} if (f) { s=t;t=[];l=s.length; }\n          }\n          var ret=(sign==0)?'0':((sign==1)?'':'-')+(t.length?'e'+t.join(''):'1'); return (brm&&brm[ret])||(brm&&brm['-'+ret]&&'-'+brm['-'+ret])||ret;\n        },\n        brm=(x=>{ var ret={}; for (var i in basis) ret[basis[i]=='1'?'1':simplify(basis[i],p,q,r)] = basis[i]; return ret; })(basis);\n        \n  // As an alternative to the string fiddling, one can also bit-fiddle. In this case the basisvectors are represented by integers with 1 bit per generator set.\n    var simplify_bits = (A,B,p2)=>{ var n=p2||(p+q+r),t=0,ab=A&B,res=A^B; if (ab&((1<<r)-1)) return [0,0]; while (n--) t^=(A=A>>1); t&=B; t^=ab>>(p+r); t^=t>>16; t^=t>>8; t^=t>>4; return [1-2*(27030>>(t&15)&1),res]; },\n        bc = (v)=>{ v=v-((v>>1)& 0x55555555); v=(v&0x33333333)+((v>>2)&0x33333333); c=((v+(v>>4)&0xF0F0F0F)*0x1010101)>>24; return c }; \n  \n  if (!options.graded && tot <= 6 || options.graded===false || options.Cayley) {\n  // Faster and degenerate-metric-resistant dualization. (a remapping table that maps items into their duals).         \n    var drm=basis.map((a,i)=>{ return {a:a,i:i} })\n                 .sort((a,b)=>a.a.length>b.a.length?1:a.a.length<b.a.length?-1:(+a.a.slice(1).split('').sort().join(''))-(+b.a.slice(1).split('').sort().join('')) )\n                 .map(x=>x.i).reverse(),\n        drms=drm.map((x,i)=>(x==0||i==0)?1:simplify(basis[x]+basis[i])[0]=='-'?-1:1);\n \n  /// Store the full metric (also for bivectors etc ..)         \n    var metric = basis.map((x,xi)=>simplify(x+x,p,q,r)|0);\n    \n  /// Generate multiplication tables for the outer and geometric products.  \n    var mulTable   = options.Cayley||basis.map(x=>basis.map(y=>(x==1)?y:(y==1)?x:simplify(x+y,p,q,r)));\n\n  /// Convert Cayley table to product matrices. The outer product selects the strict sum of the GP (but without metric), the inner product\n  /// is the left contraction.           \n    var gp=basis.map(x=>basis.map(x=>'0')), cp=gp.map(x=>gp.map(x=>'0')), cps=gp.map(x=>gp.map(x=>'0')), op=gp.map(x=>gp.map(x=>'0')), gpo={};          // Storage for our product tables.\n    basis.forEach((x,xi)=>basis.forEach((y,yi)=>{ var n = mulTable[xi][yi].replace(/^-/,''); if (!gpo[n]) gpo[n]=[]; gpo[n].push([xi,yi]); }));\n    basis.forEach((o,oi)=>{\n      gpo[o].forEach(([xi,yi])=>op[oi][xi]=(grades[oi]==grades[xi]+grades[yi])?((mulTable[xi][yi]=='0')?'0':((mulTable[xi][yi][0]!='-')?'':'-')+'b['+yi+']*this['+xi+']'):'0');\n      gpo[o].forEach(([xi,yi])=>{\n        gp[oi][xi] =((gp[oi][xi]=='0')?'':gp[oi][xi]+'+')   + ((mulTable[xi][yi]=='0')?'0':((mulTable[xi][yi][0]!='-')?'':'-')+'b['+yi+']*this['+xi+']');\n        cp[oi][xi] =((cp[oi][xi]=='0')?'':cp[oi][xi]+'+')   + ((grades[oi]==grades[yi]-grades[xi])?gp[oi][xi]:'0'); \n        cps[oi][xi]=((cps[oi][xi]=='0')?'':cps[oi][xi]+'+') + ((grades[oi]==Math.abs(grades[yi]-grades[xi]))?gp[oi][xi]:'0'); \n      });\n    });\n    \n  /// Flat Algebra Multivector Base Class.\n    var generator = class MultiVector extends (options.baseType||Float32Array) {\n    /// constructor - create a floating point array with the correct number of coefficients.\n      constructor(a) { super(a||basis.length); return this; }\n \n    /// grade selection - return a only the part of the input with the specified grade.  \n      Grade(grade,res) { res=res||new this.constructor(); for (var i=0,l=res.length; i<l; i++) if (grades[i]==grade) res[i]=this[i]; else res[i]=0; return res; }\n      Even(res) { res=res||new this.constructor(); for (var i=0,l=res.length; i<l; i++) if (grades[i]%2==0) res[i]=this[i]; else res[i]=0; return res; }\n      \n    /// grade creation - convert array with just one grade to full multivector.\n      nVector(grade,...args) { this.set(args,grade_start[grade]); return this; }\n    \n    /// Fill in coordinates (accepts sequence of index,value as arguments)\n      Coeff() { for (var i=0,l=arguments.length; i<l; i+=2) this[arguments[i]]=arguments[i+1]; return this; }\n      \n    /// Negates specific grades (passed in as args)\n      Map(res, ...a) { for (var i=0, l=res.length; i<l; i++) res[i] = (~a.indexOf(grades[i]))?-this[i]:this[i]; return res; }\n\n    /// Returns the vector grade only.\n      get Vector ()    { return this.slice(grade_start[1],grade_start[2]); };\n\n      toString() { var res=[]; for (var i=0; i<basis.length; i++) if (Math.abs(this[i])>1e-10) res.push(((this[i]==1)&&i?'':((this[i]==-1)&&i)?'-':(this[i].toFixed(10)*1))+(i==0?'':tot==1&&q==1?'i':basis[i].replace('e','e_'))); return res.join('+').replace(/\\+-/g,'-')||'0'; }\n    }  \n    \n  /// Convert symbolic matrices to code. (skipping zero's on dot and wedge matrices).\n  /// These all do straightforward string fiddling. If the 'mix' option is set they reference basis components using e.g. '.e1' instead of eg '[3]' .. so that\n  /// it will work for elements of subalgebras etc.\n    generator.prototype.Add   = new Function('b,res','res=res||new this.constructor();\\n'+basis.map((x,xi)=>'res['+xi+']=b['+xi+']+this['+xi+']').join(';\\n').replace(/(b|this)\\[(.*?)\\]/g,(a,b,c)=>options.mix?'('+b+'.'+(c|0?basis[c]:'s')+'||0)':a)+';\\nreturn res')\n    generator.prototype.Scale = new Function('b,res','res=res||new this.constructor();\\n'+basis.map((x,xi)=>'res['+xi+']=b*this['+xi+']').join(';\\n').replace(/(b|this)\\[(.*?)\\]/g,(a,b,c)=>options.mix?'('+b+'.'+(c|0?basis[c]:'s')+'||0)':a)+';\\nreturn res')\n    generator.prototype.Sub   = new Function('b,res','res=res||new this.constructor();\\n'+basis.map((x,xi)=>'res['+xi+']=this['+xi+']-b['+xi+']').join(';\\n').replace(/(b|this)\\[(.*?)\\]/g,(a,b,c)=>options.mix?'('+b+'.'+(c|0?basis[c]:'s')+'||0)':a)+';\\nreturn res')\n    generator.prototype.Mul   = new Function('b,res','res=res||new this.constructor();\\n'+gp.map((r,ri)=>'res['+ri+']='+r.join('+').replace(/\\+\\-/g,'-').replace(/(\\w*?)\\[(.*?)\\]/g,(a,b,c)=>options.mix?'('+b+'.'+(c|0?basis[c]:'s')+'||0)':a).replace(/\\+0/g,'')+';').join('\\n')+'\\nreturn res;');\n    generator.prototype.LDot  = new Function('b,res','res=res||new this.constructor();\\n'+cp.map((r,ri)=>'res['+ri+']='+r.join('+').replace(/\\+\\-/g,'-').replace(/\\+0/g,'').replace(/(\\w*?)\\[(.*?)\\]/g,(a,b,c)=>options.mix?'('+b+'.'+(c|0?basis[c]:'s')+'||0)':a)+';').join('\\n')+'\\nreturn res;');\n    generator.prototype.Dot   = new Function('b,res','res=res||new this.constructor();\\n'+cps.map((r,ri)=>'res['+ri+']='+r.join('+').replace(/\\+\\-/g,'-').replace(/\\+0/g,'').replace(/(\\w*?)\\[(.*?)\\]/g,(a,b,c)=>options.mix?'('+b+'.'+(c|0?basis[c]:'s')+'||0)':a)+';').join('\\n')+'\\nreturn res;');\n    generator.prototype.Wedge = new Function('b,res','res=res||new this.constructor();\\n'+op.map((r,ri)=>'res['+ri+']='+r.join('+').replace(/\\+\\-/g,'-').replace(/\\+0/g,'').replace(/(\\w*?)\\[(.*?)\\]/g,(a,b,c)=>options.mix?'('+b+'.'+(c|0?basis[c]:'s')+'||0)':a)+';').join('\\n')+'\\nreturn res;');\n    generator.prototype.Vee   = new Function('b,res','res=res||new this.constructor();\\n'+op.map((r,ri)=>'res['+drm[ri]+']='+r.map(x=>x.replace(/\\[(.*?)\\]/g,function(a,b){return '['+(drm[b|0])+']'})).join('+').replace(/\\+\\-/g,'-').replace(/\\+0/g,'').replace(/(\\w*?)\\[(.*?)\\]/g,(a,b,c)=>options.mix?'('+b+'.'+(c|0?basis[c]:'s')+'||0)':a)+';').join('\\n')+'\\nreturn res;');\n\n  /// Add getter and setters for the basis vectors/bivectors etc .. \n    basis.forEach((b,i)=>{generator.prototype.__defineGetter__(i?b:'s',function(){ return this[i] }); }); \n    basis.forEach((b,i)=>{generator.prototype.__defineSetter__(i?b:'s',function(x){ this[i]=x; }); });\n    \n  /// Reversion, Involutions, Conjugation for any number of grades, component acces shortcuts.\n    generator.prototype.__defineGetter__('Negative', function(){ var res = new this.constructor(); for (var i=0; i<this.length; i++) res[i]= -this[i]; return res; });\n    generator.prototype.__defineGetter__('Reverse',  function(){ var res = new this.constructor(); for (var i=0; i<this.length; i++) res[i]= this[i]*[1,1,-1,-1][grades[i]%4]; return res; });\n    generator.prototype.__defineGetter__('Involute', function(){ var res = new this.constructor(); for (var i=0; i<this.length; i++) res[i]= this[i]*[1,-1,1,-1][grades[i]%4]; return res; });\n    generator.prototype.__defineGetter__('Conjugate',function(){ var res = new this.constructor(); for (var i=0; i<this.length; i++) res[i]= this[i]*[1,-1,-1,1][grades[i]%4]; return res; });\n  \n  /// The Dual, Length, non-metric length and normalized getters.  \n    generator.prototype.__defineGetter__('Dual',function(){ if (r) return this.map((x,i,a)=>a[drm[i]]*drms[i]); var res = new this.constructor(); res[res.length-1]=1; return res.Mul(this); });\n    generator.prototype.__defineGetter__('Length',  function(){  return Math.sqrt(Math.abs(this.Mul(this.Conjugate).s)); }); \n    generator.prototype.__defineGetter__('VLength',  function(){ var res = 0; for (var i=0; i<this.length; i++) res += this[i]*this[i]; return Math.sqrt(res); });\n    generator.prototype.__defineGetter__('Normalized', function(){ var res = new this.constructor(),l=this.Length; if (!l) return this; l=1/l; for (var i=0; i<this.length; i++) res[i]=this[i]*l; return res; });\n  \n  /// Graded generator for high-dimensional algebras.\n  } else {\n  \n  /// extra graded lookups.\n    var basisg = grade_start.slice(0,grade_start.length-1).map((x,i)=>basis.slice(x,grade_start[i+1]));\n    var counts = grade_start.map((x,i,a)=>i==a.length-1?0:a[i+1]-x).slice(0,tot+1);\n    var basis_bits = basis.map(x=>x=='1'?0:x.slice(1).match(tot>9?/\\d\\d/g:/\\d/g).reduce((a,b)=>a+(1<<(b-low)),0)),\n        bits_basis = []; basis_bits.forEach((b,i)=>bits_basis[b]=i);    \n    var metric = basisg.map((x,xi)=>x.map((y,yi)=>simplify_bits(basis_bits[grade_start[xi]+yi],basis_bits[grade_start[xi]+yi])[0]));\n    var drms   = basisg.map((x,xi)=>x.map((y,yi)=>simplify_bits(basis_bits[grade_start[xi]+yi],(~basis_bits[grade_start[xi]+yi])&((1<<tot)-1))[0]));\n    \n  /// Flat Algebra Multivector Base Class.\n    var generator = class MultiVector extends Array {\n    /// constructor - create a floating point array with the correct number of coefficients.\n      constructor(a) { super(a||tot); return this; }\n \n    /// grade selection - return a only the part of the input with the specified grade.  \n      Grade(grade,res) { res=new this.constructor(); res[grade] = this[grade]; return res; }\n      \n    /// grade creation - convert array with just one grade to full multivector.\n      nVector(grade,...args) { this[grade]=args; return this; }\n    \n    /// Fill in coordinates (accepts sequence of index,value as arguments)\n      Coeff() { \n                for (var i=0,l=arguments.length; i<l; i+=2) { \n                 var gi = grades[arguments[i]];\n                 if (this[gi]==undefined) this[gi]=[];\n                 this[gi][arguments[i]-grade_start[gi]]=arguments[i+1]; \n                } \n                return this; \n              }\n      \n    /// Negates specific grades (passed in as args)\n      Map(res, ...a) {  /* tbc */ }\n\n    /// Returns the vector grade only.\n      get Vector ()    { return this[1] };\n      \n    /// multivector addition, subtraction and scalar multiplication.\n      Add(b,r) {\n        r=r||new this.constructor();\n        for (var i=0,l=Math.max(this.length,b.length);i<l;i++) \n          if (!this[i] || !b[i]) r[i] = (!this[i]) ? b[i]:this[i];\n          else { if (r[i]==undefined) r[i]=[]; for(var j=0,m=Math.max(this[i].length,b[i].length);j<m;j++) \n          {\n            if (typeof this[i][j]==\"string\" || typeof r[i][j]==\"string\" || typeof b[i][j]==\"string\") {\n             if (!this[i][j]) r[i][j] = \"\"+b[i][j];\n             else if (!b[i][j]) r[i][j] = \"\"+this[i][j];\n             else r[i][j]=\"(\"+(this[i][j]||\"0\")+(b[i][j][0]==\"-\"?\"\":\"+\")+(b[i][j]||\"0\")+\")\"; \n            } else r[i][j]=(this[i][j]||0)+(b[i][j]||0); \n          }}\n        return r;\n      }\n      Sub(b,r) {\n        r=r||new this.constructor();\n        for (var i=0,l=Math.max(this.length,b.length);i<l;i++) \n          if (!this[i] || !b[i]) r[i] = (!this[i]) ? (b[i]?b[i].map(x=>(typeof x==\"string\")?\"-\"+x:-x):undefined):this[i];\n          else { if (r[i]==undefined) r[i]=[]; for(var j=0,m=Math.max(this[i].length,b[i].length);j<m;j++) \n            if (typeof this[i][j]==\"string\" || typeof r[i][j]==\"string\" || typeof b[i][j]==\"string\") r[i][j]=\"(\"+(this[i][j]||\"0\")+\"-\"+(b[i][j]||\"0\")+\")\"; \n            else r[i][j]=(this[i][j]||0)-(b[i][j]||0); \n          }    \n        return r;\n      }\n      Scale(s) { return this.map(x=>x&&x.map(y=>typeof y==\"string\"?y+\"*\"+s:y*s)); }      \n      \n    // geometric product.\n      Mul(b,r) {\n        r=r||new this.constructor(); var gotstring=false;\n        for (var i=0,x,gsx; gsx=grade_start[i],x=this[i],i<this.length; i++) if (x) for (var j=0,y,gsy;gsy=grade_start[j],y=b[j],j<b.length; j++) if (y) for (var a=0; a<x.length; a++) if (x[a]) for (var bb=0; bb<y.length; bb++) if (y[bb]) {\n          if (i==j && a==bb) { r[0] = r[0]||(typeof x[0]==\"string\" || typeof y[bb]==\"string\"?[\"\"]:[0]); \n            if (typeof x[a]==\"string\" || typeof r[0][0]==\"string\" || typeof y[bb]==\"string\") {\n            r[0][0] = (r[0][0]?(r[0][0]+(x[a][0]==\"-\"?\"\":\"+\")):\"\")+ x[a]+\"*\"+y[bb]+(metric[i][a]!=1?\"*\"+metric[i][a]:\"\");  gotstring=true;\n            } else r[0][0] += x[a]*y[bb]*metric[i][a]; \n          } else { \n             var rn=simplify_bits(basis_bits[gsx+a],basis_bits[gsy+bb]), g=bc(rn[1]), e=bits_basis[rn[1]]-grade_start[g]; \n             if (!r[g])r[g]=[]; \n               if (typeof r[g][e]==\"string\"||typeof x[a]==\"string\"||typeof y[bb]==\"string\") { \n                 r[g][e] = (r[g][e]?r[g][e]+\"+\":\"\") + (rn[0]!=1?rn[0]+\"*\":\"\")+ x[a]+(y[bb]!=1?\"*\"+y[bb]:\"\"); gotstring=true;\n               } else r[g][e] = (r[g][e]||0) + rn[0]*x[a]*y[bb];\n          }  \n        }\n        if (gotstring) return r.map(g=>g.map(e=>e&&'('+e+')'))\n        return r;\n      }    \n    // outer product.     \n      Wedge(b,r) {\n        r=r||new this.constructor();\n        for (var i=0,x,gsx; gsx=grade_start[i],x=this[i],i<this.length; i++) if (x) for (var j=0,y,gsy;gsy=grade_start[j],y=b[j],j<b.length; j++) if (y) for (var a=0; a<x.length; a++) if (x[a]) for (var bb=0; bb<y.length; bb++) if (y[bb]) {\n          if (i!=j || a!=bb) { \n             var n1=basis_bits[gsx+a], n2=basis_bits[gsy+bb], rn=simplify_bits(n1,n2,tot), g=bc(rn[1]), e=bits_basis[rn[1]]-grade_start[g]; \n             if (g == i+j) { if (!r[g]) r[g]=[]; r[g][e] = (r[g][e]||0) + rn[0]*x[a]*y[bb]; }\n          }  \n        }\n        return r;\n      }    \n    // outer product glsl output.     \n      OPNS_GLSL(b,point_source) {\n        var r='',count=0,curg;\n        for (var i=0,x,gsx; gsx=grade_start[i],x=this[i],i<this.length; i++) if (x) for (var j=0,y,gsy;gsy=grade_start[j],y=b[j],j<b.length; j++) if (y) for (var a=0; a<counts[i]; a++) for (var bb=0; bb<counts[j]; bb++) {\n          if (i!=j || a!=bb) { \n             var n1=basis_bits[gsx+a], n2=basis_bits[gsy+bb], rn=simplify_bits(n1,n2,tot), g=bc(rn[1]), e=bits_basis[rn[1]]-grade_start[g]; \n             if (g == i+j) { curg=g; r += `res[${e}]${rn[0]=='1'?\"+=\":\"-=\"}(${point_source[a]})*b[${bb}]; //${count++}\\n`;  }\n          }  \n        }\n        r=r.split('\\n').filter(x=>x).sort((a,b)=>((a.match(/\\d+/)[0]|0)-(b.match(/\\d+/)[0]|0))||((a.match(/\\d+$/)[0]|0)-(b.match(/\\d+$/)[0]|0))).map(x=>x.replace(/\\/\\/\\d+$/,''));\n        var r2 = 'float sum=0.0; float res=0.0;\\n', g=0;\n        r.forEach(x=>{\n          var cg = x.match(/\\d+/)[0]|0;\n          if (cg != g) r2 += \"sum \"+(((metric[curg][g]==-1))?\"-=\":\"+=\")+\" res*res;\\nres = 0.0;\\n\";\n          r2 += x.replace(/\\[\\d+\\]/,'') + '\\n';\n          g=cg;\n        });\n        r2+= \"sum \"+((metric[curg][g]==-1)?\"-=\":\"+=\")+\" res*res;\\n\";\n        return r2;\n      }    \n    // Left contraction.\n      LDot(b,r) {\n        r=r||new this.constructor();\n        for (var i=0,x,gsx; gsx=grade_start[i],x=this[i],i<this.length; i++) if (x) for (var j=0,y,gsy;gsy=grade_start[j],y=b[j],j<b.length; j++) if (y) for (var a=0; a<x.length; a++) if (x[a]) for (var bb=0; bb<y.length; bb++) if (y[bb]) {\n          if (i==j && a==bb) { r[0] = r[0]||[0]; r[0][0] += x[a]*y[bb]*metric[i][a]; } \n          else { \n             var rn=simplify_bits(basis_bits[gsx+a],basis_bits[gsy+bb]), g=bc(rn[1]), e=bits_basis[rn[1]]-grade_start[g]; \n             if (g == j-i) { if (!r[g])r[g]=[]; r[g][e] = (r[g][e]||0) + rn[0]*x[a]*y[bb]; }\n          }  \n        }\n        return r;\n      }    \n    // Symmetric contraction.\n      Dot(b,r) {\n        r=r||new this.constructor();\n        for (var i=0,x,gsx; gsx=grade_start[i],x=this[i],i<this.length; i++) if (x) for (var j=0,y,gsy;gsy=grade_start[j],y=b[j],j<b.length; j++) if (y) for (var a=0; a<x.length; a++) if (x[a]) for (var bb=0; bb<y.length; bb++) if (y[bb]) {\n          if (i==j && a==bb) { r[0] = r[0]||[0]; r[0][0] += x[a]*y[bb]*metric[i][a]; } \n          else { \n             var rn=simplify_bits(basis_bits[gsx+a],basis_bits[gsy+bb]), g=bc(rn[1]), e=bits_basis[rn[1]]-grade_start[g]; \n             if (g == Math.abs(j-i)) { if (!r[g])r[g]=[]; r[g][e] = (r[g][e]||0) + rn[0]*x[a]*y[bb]; }\n          }  \n        }\n        return r;\n      }    \n    // Should be optimized..  \n      Vee(b,r)          { return (this.Dual.Wedge(b.Dual)).Dual; }\n    // Output, lengths, involutions, normalized, dual.  \n      toString()        { return [...this].map((g,gi)=>g&&g.map((c,ci)=>!c?undefined:c+basisg[gi][ci]).filter(x=>x).join('+')).filter(x=>x).join('+').replace(/\\+\\-/g,'-'); }  \n      get s ()          { if (this[0]) return this[0][0]||0; return 0; }\n      get Length ()     { var res=0; this.forEach((g,gi)=>g&&g.forEach((e,ei)=>res+=(e||0)**2*metric[gi][ei])); return Math.abs(res)**.5; }\n      get VLength ()    { var res=0; this.forEach((g,gi)=>g&&g.forEach((e,ei)=>res+=(e||0)**2)); return Math.abs(res)**.5; }\n      get Reverse ()    { var r=new this.constructor(); this.forEach((x,gi)=>x&&x.forEach((e,ei)=>{if(!r[gi])r[gi]=[]; r[gi][ei] = this[gi][ei]*[1,1,-1,-1][gi%4]; })); return r; }\n      get Involute ()   { var r=new this.constructor(); this.forEach((x,gi)=>x&&x.forEach((e,ei)=>{if(!r[gi])r[gi]=[]; r[gi][ei] = this[gi][ei]*[1,-1,1,-1][gi%4]; })); return r; }\n      get Conjugate ()  { var r=new this.constructor(); this.forEach((x,gi)=>x&&x.forEach((e,ei)=>{if(!r[gi])r[gi]=[]; r[gi][ei] = this[gi][ei]*[1,-1,-1,1][gi%4]; })); return r; }\n      get Dual()        { var r=new this.constructor(); this.forEach((g,gi)=>{ if (!g) return; r[tot-gi]=[]; g.forEach((e,ei)=>r[tot-gi][counts[gi]-1-ei]=drms[gi][ei]*e); }); return r; }\n      get Normalized () { return this.Scale(1/this.Length); }\n    }  \n\n    \n   // This generator is UNDER DEVELOPMENT - I'm publishing it so I can test on observable.\n  }\n  \n  // Generate a new class for our algebra. It extends the javascript typed arrays (default float32 but can be specified in options).\n    var res = class Element extends generator {\n    \n    // constructor - create a floating point array with the correct number of coefficients.\n      constructor(a) { super(a); return this; }\n      \n    // Grade selection. (implemented by parent class).   \n      Grade(grade,res) { res=res||new Element(); return super.Grade(grade,res); }\n      \n    // Right and Left divide - Defined on the elements, shortcuts to multiplying with the inverse.  \n      Div  (b,res) { return this.Mul(b.Inverse,res); }\n      LDiv (b,res) { return b.Inverse.Mul(this,res); }\n    \n    // Taylor exp - I will replace this with something smarter for elements of the even subalgebra's and other pure blades.  \n      Exp  ()      { \n        if (r==1 && tot<=4 && this[0]==0) { \n          var sq = (tot==4)?-(this[8]**2+this[9]**2+this[10]**2):this.Mul(this).s;  if (sq==0) { var res = this.slice();res[0]+=1; return res; }\n          var l = Math.sqrt(Math.abs(sq)); if (sq<0)  { var res = this.Scale( Math.sin(l)/l ); res[0]=Math.cos(l); return res; }\n          var res = this.Scale( Math.sinh(l)/l ); res[0]=Math.cosh(l); return res;\n        }\n        var res = Element.Scalar(1), y=1, M= new Element(this), N=new Element(this); for (var x=1; x<25; x++) { res=res.Add(M.Mul(Element.Scalar(1/y))); M=M.Mul(N); y=y*(x+1); }; return res; \n      }\n      \n    // Helper for efficient inverses. (custom involutions - negates grades in arguments). \n      Map () { var res=new Element(); return super.Map(res,...arguments); }\n      \n    // Factories - Make it easy to generate vectors, bivectors, etc when using the functional API. None of the examples use this but\n    // users that have used other GA libraries will expect these calls. The Coeff() is used internally when translating algebraic literals.\n      static Element()   { return new Element([...arguments]); };\n      static Coeff()     { return (new Element()).Coeff(...arguments); }\n      static Scalar(x)   { return (new Element()).Coeff(0,x); }\n      static Vector()    { return (new Element()).nVector(1,...arguments); }\n      static Bivector()  { return (new Element()).nVector(2,...arguments); }\n      static Trivector() { return (new Element()).nVector(3,...arguments); }\n      static nVector(n)  { return (new Element()).nVector(...arguments); }\n    \n    // Static operators. The parser will always translate operators to these static calls so that scalars, vectors, matrices and other non-multivectors can also be handled.\n    // The static operators typically handle functions and matrices, calling through to element methods for multivectors. They are intended to be flexible and allow as many\n    // types of arguments as possible. If performance is a consideration, one should use the generated element methods instead. (which only accept multivector arguments)\n      static toEl(x)        { if (x instanceof Function) x=x(); if (!(x instanceof Element)) x=Element.Scalar(x); return x; }\n    \n    // Addition and subtraction. Subtraction with only one parameter is negation.   \n      static Add(a,b,res)   { \n      // Resolve expressions passed in.\n        while(a.call)a=a(); while(b.call)b=b(); if (a.Add && b.Add) return a.Add(b,res);  \n      // If either is a string, the result is a string.  \n        if ((typeof a=='string')||(typeof b=='string')) return a.toString()+b.toString(); \n      // If only one is an array, add the other element to each of the elements.   \n        if ((a instanceof Array)^(b instanceof Array)) return (a instanceof Array)?a.map(x=>Element.Add(x,b)):b.map(x=>Element.Add(a,x)); \n      // If both are equal length arrays, add elements one-by-one  \n        if ((a instanceof Array)&&(b instanceof Array)&&a.length==b.length) return a.map((x,xi)=>Element.Add(x,b[xi])); \n      // If they're both not elements let javascript resolve it.  \n        if (!(a instanceof Element || b instanceof Element)) return a+b; \n      // Here we're left with scalars and multivectors, call through to generated code.   \n        a=Element.toEl(a); b=Element.toEl(b); return a.Add(b,res);\n      }\n      \n      static Sub(a,b,res)   {  \n      // Resolve expressions passed in.\n        while(a.call)a=a(); while(b&&b.call) b=b(); if (a.Sub && b && b.Sub) return a.Sub(b,res);\n      // If only one is an array, add the other element to each of the elements.   \n        if (b&&((a instanceof Array)^(b instanceof Array))) return (a instanceof Array)?a.map(x=>Element.Sub(x,b)):b.map(x=>Element.Sub(a,x)); \n      // If both are equal length arrays, add elements one-by-one  \n        if (b&&(a instanceof Array)&&(b instanceof Array)&&a.length==b.length) return a.map((x,xi)=>Element.Sub(x,b[xi])); \n      // Negation  \n        if (arguments.length==1) return Element.Mul(a,-1); \n      // If none are elements here, let js do it.  \n        if (!(a instanceof Element || b instanceof Element)) return a-b; \n      // Here we're left with scalars and multivectors, call through to generated code.   \n        a=Element.toEl(a); b=Element.toEl(b); return a.Sub(b,res);\n      }\n    \n    // The geometric product. (or matrix*matrix, matrix*vector, vector*vector product if called with 1D and 2D arrays)\n      static Mul(a,b,res)   {\n      // Resolve expressions  \n        while(a.call&&!a.length)a=a(); while(b.call&&!b.length)b=b(); if (a.Mul && b.Mul) return a.Mul(b,res);\n      // still functions -> experimental curry style (dont use this.)\n        if (a.call && b.call) return (ai,bi)=>Element.Mul(a(ai),b(bi));   \n      // scalar mul.\n        if (Number.isFinite(a) && b.Scale) return b.Scale(a); else if (Number.isFinite(b) && a.Scale) return a.Scale(b);  \n      // Handle matrices and vectors.  \n        if ((a instanceof Array)&&(b instanceof Array)) { \n        // vector times vector performs a dot product. (which internally uses the GP on each component)\n          if((!(a[0] instanceof Array) || (a[0] instanceof Element)) &&(!(b[0] instanceof Array) || (b[0] instanceof Element))) { var r=tot?Element.Scalar(0):0; a.forEach((x,i)=>r=Element.Add(r,Element.Mul(x,b[i]),r)); return r; } \n        // Array times vector  \n          if(!(b[0] instanceof Array)) return a.map((x,i)=>Element.Mul(a[i],b)); \n        // Array times Array  \n          var r=a.map((x,i)=>b[0].map((y,j)=>{ var r=tot?Element.Scalar(0):0; x.forEach((xa,k)=>r=Element.Add(r,Element.Mul(xa,b[k][j]))); return r; })); \n        // Return resulting array or scalar if 1 by 1.   \n          if (r.length==1 && r[0].length==1) return r[0][0]; else return r;  \n        } \n      // Only one is an array multiply each of its elements with the other.  \n        if ((a instanceof Array)^(b instanceof Array)) return (a instanceof Array)?a.map(x=>Element.Mul(x,b)):b.map(x=>Element.Mul(a,x)); \n      // Try js multiplication, else call through to geometric product.  \n        var r=a*b; if (!isNaN(r)) return r; \n        a=Element.toEl(a); b=Element.toEl(b); return a.Mul(b,res);\n      }  \n      \n    // The inner product. (default is left contraction).  \n      static LDot(a,b,res)   {  \n      // Expressions\n        while(a.call)a=a(); while(b.call)b=b(); //if (a.LDot) return a.LDot(b,res);\n      // Map elements in array  \n        if (b instanceof Array && !(a instanceof Array)) return b.map(x=>Element.LDot(a,x)); \n        if (a instanceof Array && !(b instanceof Array)) return a.map(x=>Element.LDot(x,b)); \n      // js if numbers, else contraction product.  \n        if (!(a instanceof Element || b instanceof Element)) return a*b; \n        a=Element.toEl(a);b=Element.toEl(b); return a.LDot(b,res); \n      }  \n \n    // The symmetric inner product. (default is left contraction).  \n      static Dot(a,b,res)   {  \n      // Expressions\n        while(a.call)a=a(); while(b.call)b=b(); //if (a.LDot) return a.LDot(b,res);\n      // js if numbers, else contraction product.  \n        if (!(a instanceof Element || b instanceof Element)) return a|b; \n        a=Element.toEl(a);b=Element.toEl(b); return a.Dot(b,res); \n      }  \n      \n    // The outer product. (Grassman product - no use of metric)  \n      static Wedge(a,b,res) {  \n      // Expressions\n        while(a.call)a=a(); while(b.call)b=b(); if (a.Wedge) return a.Wedge(Element.toEl(b),res); \n      // The outer product of two vectors is a matrix .. internally Mul not Wedge !  \n        if (a instanceof Array && b instanceof Array) return a.map(xa=>b.map(xb=>Element.Mul(xa,xb)));\n      // js, else generated wedge product.\n        if (!(a instanceof Element || b instanceof Element)) return a*b; \n        a=Element.toEl(a);b=Element.toEl(b); return a.Wedge(b,res); \n      }  \n      \n    // The regressive product. (Dual of the outer product of the duals). \n      static Vee(a,b,res) {  \n      // Expressions\n        while(a.call)a=a(); while(b.call)b=b(); if (a.Vee) return a.Vee(Element.toEl(b),res);\n      // js, else generated vee product. (shortcut for dual of wedge of duals)\n        if (!(a instanceof Element || b instanceof Element)) return 0; \n        a=Element.toEl(a);b=Element.toEl(b); return a.Vee(b,res); \n      }  \n     \n    // The sandwich product. Provided for convenience (>>> operator)  \n      static sw(a,b) {  \n      // Expressions\n        while(a.call)a=a(); while(b.call)b=b(); if (a.sw) return a.sw(b);\n      // Map elements in array  \n        if (b instanceof Array) return b.map(x=>Element.sw(a,x)); \n      // Call through. no specific generated code for it so just perform the muls.  \n        a=Element.toEl(a); b=Element.toEl(b); return a.Mul(b).Mul(a.Conjugate); \n      }\n\n    // Division - scalars or cal through to element method.\n      static Div(a,b,res) {  \n      // Expressions\n        while(a.call)a=a(); while(b.call)b=b();\n      // js or call through to element divide.  \n        if (!(a instanceof Element || b instanceof Element)) return a/b; \n        a=Element.toEl(a);\n        if (Number.isFinite(b)) { return a.Scale(1/b,res); }\n        b=Element.toEl(b); return a.Div(b,res); \n      }  \n      \n    // Pow - needs obvious extensions for natural powers. (exponentiation by squaring)  \n      static Pow(a,b,res) {  \n      // Expressions\n        while(a.call)a=a(); while(b.call)b=b(); if (a.Pow) return a.Pow(b,res); \n      // Exponentiation.  \n        if (a==Math.E && b.Exp) return b.Exp(); \n      // Squaring  \n        if (b===2) return this.Mul(a,a,res);\n      // No elements, call through to js  \n        if (!(a instanceof Element || b instanceof Element)) return a**b; \n      // Inverse  \n        if (b===-1) return a.Inverse;  \n      // Call through to element pow.  \n        a=Element.toEl(a); return a.Pow(b); \n      }  \n      \n    // Handles scalars and calls through to element method.  \n      static exp(a) {  \n      // Expressions.\n        while(a.call)a=a(); \n      // If it has an exp callthrough, use it, else call through to math.  \n        if (a.Exp) return a.Exp(); \n        return Math.exp(a); \n      }\n      \n    // Dual, Involute, Reverse, Conjugate, Normalize and length, all direct call through. Conjugate handles matrices.\n      static Dual(a)      { return Element.toEl(a).Dual; }; \n      static Involute(a)  { return Element.toEl(a).Involute; }; \n      static Reverse(a)   { return Element.toEl(a).Reverse; }; \n      static Conjugate(a) { if (a.Conjugate) return a.Conjugate; if (a instanceof Array) return a[0].map((c,ci)=>a.map((r,ri)=>Element.Conjugate(a[ri][ci]))); return Element.toEl(a).Conjugate; }\n      static Normalize(a) { return Element.toEl(a).Normalized; }; \n      static Length(a)    { return Element.toEl(a).Length };\n      \n    // Comparison operators always use length. Handle expressions, then js or length comparison  \n      static eq(a,b)  { if (!(a instanceof Element)||!(b instanceof Element)) return a==b; while(a.call)a=a(); while(b.call)b=b(); for (var i=0; i<a.length; i++) if (a[i]!=b[i]) return false; return true; }\n      static neq(a,b) { if (!(a instanceof Element)||!(b instanceof Element)) return a!=b; while(a.call)a=a(); while(b.call)b=b(); for (var i=0; i<a.length; i++) if (a[i]!=b[i]) return true; return false; }\n      static lt(a,b)  { while(a.call)a=a(); while(b.call)b=b(); return (a instanceof Element?a.Length:a)<(b instanceof Element?b.Length:b); }\n      static gt(a,b)  { while(a.call)a=a(); while(b.call)b=b(); return (a instanceof Element?a.Length:a)>(b instanceof Element?b.Length:b); }\n      static lte(a,b) { while(a.call)a=a(); while(b.call)b=b(); return (a instanceof Element?a.Length:a)<=(b instanceof Element?b.Length:b); }\n      static gte(a,b) { while(a.call)a=a(); while(b.call)b=b(); return (a instanceof Element?a.Length:a)>=(b instanceof Element?b.Length:b); }\n      \n    // Debug output and printing multivectors.  \n      static describe(x) { if (x===true) console.log(`Basis\\n${basis}\\nMetric\\n${metric.slice(1,1+tot)}\\nCayley\\n${mulTable.map(x=>(x.map(x=>('           '+x).slice(-2-tot)))).join('\\n')}\\nMatrix Form:\\n`+gp.map(x=>x.map(x=>x.match(/(-*b\\[\\d+\\])/)).map(x=>x&&((x[1].match(/-/)||' ')+String.fromCharCode(65+1*x[1].match(/\\d+/)))||' 0')).join('\\n')); return {basis:basisg||basis,metric,mulTable} }    \n      \n    // Direct sum of algebras - experimental\n      static sum(B){\n        var A = Element;\n        // Get the multiplication tabe and basis.\n        var T1 = A.describe().mulTable, T2 = B.describe().mulTable;\n        var B1 = A.describe().basis, B2 = B.describe().basis;\n        // Get the maximum index of T1, minimum of T2 and rename T2 if needed.\n        var max_T1 = B1.filter(x=>x.match(/e/)).map(x=>x.match(/\\d/g)).flat().map(x=>x|0).sort((a,b)=>b-a)[0];\n        var max_T2 = B2.filter(x=>x.match(/e/)).map(x=>x.match(/\\d/g)).flat().map(x=>x|0).sort((a,b)=>b-a)[0];\n        var min_T2 = B2.filter(x=>x.match(/e/)).map(x=>x.match(/\\d/g)).flat().map(x=>x|0).sort((a,b)=>a-b)[0];\n        // remapping ..\n        T2 = T2.map(x=>x.map(y=>y.match(/e/)?y.replace(/(\\d)/g,(x)=>(x|0)+max_T1):y.replace(\"1\",\"e\"+(1+max_T2+max_T1))));\n        B2 = B2.map((y,i)=>i==0?y.replace(\"1\",\"e\"+(1+max_T2+max_T1)):y.replace(/(\\d)/g,(x)=>(x|0)+max_T1));\n        // Build the new basis and multable..\n        var basis = [...B1,...B2];\n        var Cayley = T1.map((x,i)=>[...x,...T2[0].map(x=>\"0\")]).concat(T2.map((x,i)=>[...T1[0].map(x=>\"0\"),...x]))\n        // Build the new algebra.\n        var grades = [...B1.map(x=>x==\"1\"?0:x.length-1),...B2.map((x,i)=>i?x.length-1:0)];\n        var a = Algebra({basis,Cayley,grades,tot:Math.log2(B1.length)+Math.log2(B2.length)})\n        // And patch up ..\n        a.Scalar = function(x) {\n          var res = new a();\n          for (var i=0; i<res.length; i++) res[i] = basis[i] == Cayley[i][i] ? x:0;\n          return res;\n        }\n        return a;\n      }  \n      \n    // The graphing function supports several modes. It can render 1D functions and 2D functions on canvas, and PGA2D, PGA3D and CGA2D functions using SVG.\n    // It handles animation and interactivity.\n    //   graph(function(x))     => function of 1 parameter will be called with that parameter from -1 to 1 and graphed on a canvas. Returned values should also be in the [-1 1] range\n    //   graph(function(x,y))   => functions of 2 parameters will be called from -1 to 1 on both arguments. Returned values can be 0-1 for greyscale or an array of three RGB values.\n    //   graph(array)           => array of algebraic elements (points, lines, circles, segments, texts, colors, ..) is graphed.\n    //   graph(function=>array) => same as above, for animation scenario's this function is called each frame.\n    // An optional second parameter is an options object { width, height, animate, camera, scale, grid, canvas } \n      static graph(f,options) { \n      // Store the original input\n        if (!f) return; var origf=f; \n      // generate default options.  \n        options=options||{}; options.scale=options.scale||1; options.camera=options.camera||(tot<4?Element.Scalar(1):new Element([0.7071067690849304, 0, 0, 0, 0, 0, 0, 0, 0, 0.7071067690849304, 0, 0, 0, 0, 0, 0])); \n        var ww=options.width, hh=options.height, cvs=options.canvas, tpcam=new Element([0,0,0,0,0,0,0,0,0,0,0,-5,0,0,1,0]),tpy=this.Coeff(4,1),tp=new Element(), \n      // project 3D to 2D. This allows to render 3D and 2D PGA with the same code.    \n        project=(o)=>{ if (!o) return o; while (o.call) o=o(); return (tot==4 && (o.length==16))?(tpcam).Vee(options.camera.Mul(o).Mul(options.camera.Conjugate)).Wedge(tpy):(o.length==2**tot)?Element.sw(options.camera,o):o;};\n      // gl escape.\n        if (options.gl) return Element.graphGL(f,options); if (options.up) return Element.graphGL2(f,options);\n      // if we get an array or function without parameters, we render c2d or p2d SVG points/lines/circles/etc\n        if (!(f instanceof Function) || f.length===0) { \n        // Our current cursor, color, animation state and 2D mapping.\n          var lx,ly,lr,color,res,anim=false,to2d=(tot==3)?[0,1,2,3,4,5,6,7]:[0,7,9,10,13,12,14,15];\n        // Make sure we have an array of elements. (if its an object, convert to array with elements and names.)   \n          if (f instanceof Function) f=f(); if (!(f instanceof Array)) f=[].concat.apply([],Object.keys(f).map((k)=>typeof f[k]=='number'?[f[k]]:[f[k],k])); \n        // The build function generates the actual SVG. It will be called everytime the user interacts or the anim flag is set.  \n          function build(f,or) {\n          // Make sure we have an aray. \n            if (or && f && f instanceof Function) f=f(); \n          // Reset position and color for cursor.  \n            lx=-2;ly=-1.85;lr=0;color='#444'; \n          // Create the svg element. (master template string till end of function)  \n            var svg=new DOMParser().parseFromString(`<SVG onmousedown=\"if(evt.target==this)this.sel=undefined\" viewBox=\"-2 -${2*(hh/ww||1)} 4 ${4*(hh/ww||1)}\" style=\"width:${ww||512}px; height:${hh||512}px; background-color:#eee; -webkit-user-select:none; -moz-user-select:none; -ms-user-select:none; user-select:none\">\n            // Add a grid (option)\n            ${options.grid?[...Array(21)].map((x,xi)=>`<line x1=\"-10\" y1=\"${((xi-10)/2-(tot<4?2*options.camera.e02:0))*options.scale}\" x2=\"10\" y2=\"${((xi-10)/2-(tot<4?2*options.camera.e02:0))*options.scale}\" stroke-width=\"0.005\" stroke=\"#CCC\"/><line y1=\"-10\" x1=\"${((xi-10)/2-(tot<4?2*options.camera.e01:0))*options.scale}\" y2=\"10\" x2=\"${((xi-10)/2-(tot<4?2*options.camera.e01:0))*options.scale}\"  stroke-width=\"0.005\" stroke=\"#CCC\"/>`):''}\n            // Handle conformal 2D elements. \n            ${options.conformal?f.map&&f.map((o,oidx)=>{ \n            // Optional animation handling.\n              if((o==Element.graph && or!==false)||(oidx==0&&options.animate&&or!==false)) { anim=true; requestAnimationFrame(()=>{var r=build(origf,(!res)||(document.body.contains(res))).innerHTML; if (res) res.innerHTML=r; }); if (!options.animate) return; }\n            // Resolve expressions passed in.  \n              while (o.call) o=o();\n            // Arrays are rendered as segments or polygons. (2 or more elements)  \n              if (o instanceof Array)  { lx=ly=lr=0; o=o.map(o=>{ while(o.call)o=o(); return o; }); o.forEach((o)=>{lx+=o.e1;ly+=-o.e2});lx/=o.length;ly/=o.length; return o.length>2?`<POLYGON STYLE=\"pointer-events:none; fill:${color};opacity:0.7\" points=\"${o.map(o=>(o.e1+','+(-o.e2)+' '))}\"/>`:`<LINE style=\"pointer-events:none\" x1=${o[0].e1} y1=${-o[0].e2} x2=${o[1].e1} y2=${-o[1].e2} stroke-width=\"${options.lineWidth*0.005||0.005}\" stroke=\"${color||'#888'}\"/>`; }\n            // Strings are rendered at the current cursor position.  \n              if (typeof o =='string') { var res2=(o[0]=='_')?'':`<text x=\"${lx}\" y=\"${ly}\" font-family=\"Verdana\" font-size=\"${options.fontSize*0.1||0.1}\" style=\"pointer-events:none\" fill=\"${color||'#333'}\" transform=\"rotate(${lr},${lx},${ly})\">&nbsp;${o}&nbsp;</text>`; ly+=0.14; return res2; }\n            // Numbers change the current color.  \n              if (typeof o =='number') { color='#'+(o+(1<<25)).toString(16).slice(-6); return ''; };\n            // All other elements are rendered ..  \n              var b1=o.Grade(1).VLength>0.001,b2=o.Grade(2).VLength>0.001,b3=o.Grade(3).VLength>0.001; \n            // Points  \n              if (b1 && !b2 && !b3) { lx=o.e1; ly=-o.e2; lr=0; return res2=`<CIRCLE onmousedown=\"this.parentElement.sel=${oidx}\" cx=\"${lx}\" cy=\"${ly}\" r=\"${options.pointRadius*0.03||0.03}\" fill=\"${color||'green'}\"/>`; }\n              else if (!b1 && !b2 && b3) { var isLine=Element.Coeff(4,1,3,-1).LDot(o).Length==0; \n              // Lines.\n                if (isLine) { var loc=((Element.Coeff(4,-.5).Add(Element.Coeff(3,-.5))).LDot(o)).Div(o), att=(Element.Coeff(4,1,3,-1)).LDot(o); lx=-loc.e1; ly=loc.e2; lr=Math.atan2(att[8],att[7])/Math.PI*180; return `<LINE style=\"pointer-events:none\" x1=${lx-10} y1=${ly} x2=${lx+10} y2=${ly} stroke-width=\"${options.lineWidth*0.005||0.005}\" stroke=\"${color||'#888'}\" transform=\"rotate(${lr},${lx},${ly})\"/>`;};\n              // Circles.  \n                var loc=o.Div((Element.Coeff(4,1,3,-1)).LDot(o)); lx=-loc.e1; ly=loc.e2; var r=-o.Mul(o.Conjugate).s/(Element.Pow((Element.Coeff(4,1,3,-1)).LDot(o),2).s); r=r**0.5; return `<CIRCLE onmousedown=\"this.parentElement.sel=${oidx}\" cx=\"${lx}\" cy=\"${ly}\" r=\"${r}\" stroke-width=\"${options.lineWidth*0.005||0.005}\" fill=\"none\" stroke=\"${color||'green'}\"/>`;   \n              } else if (!b1 && b2 &&!b3) { \n              // Point Pairs.\n                lr=0; var ei=Element.Coeff(4,1,3,-1),eo=Element.Coeff(4,.5,3,.5), nix=o.Wedge(ei), sqr=o.LDot(o).s/nix.LDot(nix).s, r=Math.sqrt(Math.abs(sqr)), attitude=((ei.Wedge(eo)).LDot(nix)).Normalized.Mul(Element.Scalar(r)), pos=o.Div(nix); pos=pos.Div( pos.LDot(Element.Sub(ei))); \n                lx=pos.e1; ly=-pos.e2; if (sqr<0) return `<CIRCLE onmousedown=\"this.parentElement.sel=${oidx}\" cx=\"${lx}\" cy=\"${ly}\" r=\"${options.pointRadius*0.03||0.03}\" stroke-width=\"0.005\" fill=\"none\" stroke=\"${color||'green'}\"/>`;\n                lx=pos.e1+attitude.e1; ly=-pos.e2-attitude.e2; var res2=`<CIRCLE onmousedown=\"this.parentElement.sel=${oidx}\" cx=\"${lx}\" cy=\"${ly}\" r=\"${options.pointRadius*0.03||0.03}\" fill=\"${color||'green'}\"/>`;\n                lx=pos.e1-attitude.e1; ly=-pos.e2+attitude.e2; return res2+`<CIRCLE onmousedown=\"this.parentElement.sel=${oidx}\" cx=\"${lx}\" cy=\"${ly}\" r=\"${options.pointRadius*0.03||0.03}\" fill=\"${color||'green'}\"/>`;\n              }\n            // Handle projective 2D and 3D elements.  \n            }):f.map&&f.map((o,oidx)=>{  if((o==Element.graph && or!==false)||(oidx==0&&options.animate&&or!==false)) { anim=true; requestAnimationFrame(()=>{var r=build(origf,(!res)||(document.body.contains(res))).innerHTML; if (res) res.innerHTML=r; }); if (!options.animate) return; } while (o instanceof Function) o=o(); o=(o instanceof Array)?o.map(project):project(o); if (o===undefined) return; \n            // line segments and polygons\n              if (o instanceof Array && o.length)  { lx=ly=lr=0; o.forEach((o)=>{while (o.call) o=o(); lx+=options.scale*((drm[1]==6||drm[1]==14)?-1:1)*o[drm[2]]/o[drm[1]];ly+=options.scale*o[drm[3]]/o[drm[1]]});lx/=o.length;ly/=o.length; return o.length>2?`<POLYGON STYLE=\"pointer-events:none; fill:${color};opacity:0.7\" points=\"${o.map(o=>((drm[1]==6||drm[1]==14)?-1:1)*options.scale*o[drm[2]]/o[drm[1]]+','+options.scale*o[drm[3]]/o[drm[1]]+' ')}\"/>`:`<LINE style=\"pointer-events:none\" x1=${options.scale*((drm[1]==6||drm[1]==14)?-1:1)*o[0][drm[2]]/o[0][drm[1]]} y1=${options.scale*o[0][drm[3]]/o[0][drm[1]]} x2=${options.scale*((drm[1]==6||drm[1]==14)?-1:1)*o[1][drm[2]]/o[1][drm[1]]} y2=${options.scale*o[1][drm[3]]/o[1][drm[1]]} stroke-width=\"${options.lineWidth*0.005||0.005}\" stroke=\"${color||'#888'}\"/>`; }\n            // svg\n              if (typeof o =='string' && o[0]=='<') { return o; }  \n            // Labels  \n              if (typeof o =='string') { var res2=(o[0]=='_')?'':`<text x=\"${lx}\" y=\"${ly}\" font-family=\"Verdana\" font-size=\"${options.fontSize*0.1||0.1}\" style=\"pointer-events:none\" fill=\"${color||'#333'}\" transform=\"rotate(${lr},0,0)\">&nbsp;${o}&nbsp;</text>`; ly+=0.14; return res2; }\n            // Colors  \n              if (typeof o =='number') { color='#'+(o+(1<<25)).toString(16).slice(-6); return ''; };\n            // Points  \n              if (o[to2d[6]]**2        >0.0001) { lx=options.scale*o[drm[2]]/o[drm[1]]; if (drm[1]==6||drm[1]==14) lx*=-1; ly=options.scale*o[drm[3]]/o[drm[1]]; lr=0;  var res2=`<CIRCLE onmousedown=\"this.parentElement.sel=${oidx}\" cx=\"${lx}\" cy=\"${ly}\" r=\"${options.pointRadius*0.03||0.03}\" fill=\"${color||'green'}\"/>`; ly-=0.05; lx-=0.1; return res2; }\n            // Lines  \n              if (o[to2d[2]]**2+o[to2d[3]]**2>0.0001) { var l=Math.sqrt(o[to2d[2]]**2+o[to2d[3]]**2); o[to2d[2]]/=l; o[to2d[3]]/=l; o[to2d[1]]/=l; lx=0.5; ly=options.scale*((drm[1]==6)?-1:-1)*o[to2d[1]]; lr=-Math.atan2(o[to2d[2]],o[to2d[3]])/Math.PI*180; var res2=`<LINE style=\"pointer-events:none\" x1=-10 y1=${ly} x2=10 y2=${ly} stroke-width=\"${options.lineWidth*0.005||0.005}\" stroke=\"${color||'#888'}\" transform=\"rotate(${lr},0,0)\"/>`; ly-=0.05; return res2; }\n            // Vectors   \n              if (o[to2d[4]]**2+o[to2d[5]]**2>0.0001) { lr=0; ly+=0.05; lx+=0.1; var res2=`<LINE style=\"pointer-events:none\" x1=${lx} y1=${ly} x2=${lx-o.e02} y2=${ly+o.e01} stroke-width=\"0.005\" stroke=\"${color||'#888'}\"/>`; ly=ly+o.e01/4*3-0.05; lx=lx-o.e02/4*3; return res2; }\n            }).join()}`,'text/html').body; \n          // return the inside of the created svg element.  \n            return svg.removeChild(svg.firstChild); \n          };\n        // Create the initial svg and install the mousehandlers.  \n          res=build(f); res.value=f; res.options=options;\n          res.onmousemove=(e)=>{ if (res.sel===undefined || !e.buttons) return;var resx=res.getBoundingClientRect().width,resy=res.getBoundingClientRect().height,x=((e.clientX-res.getBoundingClientRect().left)/(resx/4||128)-2)*(resx>resy?resx/resy:1),y=((e.clientY-res.getBoundingClientRect().top)/(resy/4||128)-2)*(resy>resx?resy/resx:1);x/=options.scale;y/=options.scale; if (options.conformal) {f[res.sel][1]=x; f[res.sel][2]=-y; var l=x*x+y*y; f[res.sel][3]=0.5-l*0.5; f[res.sel][4]=0.5+l*0.5; } else {f[res.sel][drm[2]]=((drm[1]==6)?-x:x)-((tot<4)?2*options.camera.e01:0); f[res.sel][drm[3]]=y+((tot<4)?2*options.camera.e02:0); f[res.sel][drm[1]]=1;} if (!anim) res.innerHTML=build(f).innerHTML; res.dispatchEvent(new CustomEvent('input')) }; \n          return res;\n        }  \n      // 1d and 2d functions are rendered on a canvas.   \n        cvs=cvs||document.createElement('canvas'); if(ww)cvs.width=ww; if(hh)cvs.height=hh; var w=cvs.width,h=cvs.height,context=cvs.getContext('2d'), data=context.getImageData(0,0,w,h);\n      // two parameter functions .. evaluate for both and set resulting color.  \n        if (f.length==2) for (var px=0; px<w; px++) for (var py=0; py<h; py++) { var res=f(px/w*2-1, py/h*2-1); res=res.buffer?[].slice.call(res):res.slice?res:[res,res,res]; data.data.set(res.map(x=>x*255).concat([255]),py*w*4+px*4); }\n      // one parameter function.. go over x range, use result as y.   \n        else if (f.length==1) for (var px=0; px<w; px++) { var res=f(px/w*2-1); res=Math.round((res/2+0.5)*h); if (res > 0 && res < h-1) data.data.set([0,0,0,255],res*w*4+px*4); }\n        return context.putImageData(data,0,0),cvs;       \n      }\n      \n    // webGL2 Graphing function. (for OPNS/IPNS implicit 2D and 1D surfaces in 3D space).\n      static graphGL2(f,options) {\n      // Create canvas, get webGL2 context.\n        var canvas=document.createElement('canvas'); canvas.style.width=options.width||''; canvas.style.height=options.height||''; canvas.style.backgroundColor='#EEE';\n        if (options.width && options.width.match && options.width.match(/px/i)) canvas.width = parseFloat(options.width); if (options.height && options.height.match && options.height.match(/px/i)) canvas.height = parseFloat(options.height);\n        var gl=canvas.getContext('webgl2',{alpha:options.alpha||false,preserveDrawingBuffer:true,antialias:true,powerPreference:'high-performance'});\n        var gl2=!!gl; if (!gl) gl=canvas.getContext('webgl',{alpha:options.alpha||false,preserveDrawingBuffer:true,antialias:true,powerPreference:'high-performance'});\n        gl.clearColor(240/255,240/255,240/255,1.0); gl.enable(gl.DEPTH_TEST); if (!gl2) { gl.getExtension(\"EXT_frag_depth\"); gl.va = gl.getExtension('OES_vertex_array_object'); }\n        else gl.va = { createVertexArrayOES : gl.createVertexArray.bind(gl), bindVertexArrayOES : gl.bindVertexArray.bind(gl), deleteVertexArrayOES : gl.deleteVertexArray.bind(gl)  }\n      // Compile vertex and fragment shader, return program.\n        var compile=(vs,fs)=>{\n          var s=[gl.VERTEX_SHADER,gl.FRAGMENT_SHADER].map((t,i)=>{\n            var r=gl.createShader(t); gl.shaderSource(r,[vs,fs][i]); gl.compileShader(r);\n            return gl.getShaderParameter(r, gl.COMPILE_STATUS)&&r||console.error(gl.getShaderInfoLog(r));\n          });\n          var p = gl.createProgram(); gl.attachShader(p, s[0]); gl.attachShader(p, s[1]); gl.linkProgram(p);\n          gl.getProgramParameter(p, gl.LINK_STATUS)||console.error(gl.getProgramInfoLog(p));\n          return p;\n        };\n      // Create vertex array and buffers, upload vertices and optionally texture coordinates.\n        var createVA=function(vtx) {\n          var r = gl.va.createVertexArrayOES(); gl.va.bindVertexArrayOES(r);\n          var b = gl.createBuffer(); gl.bindBuffer(gl.ARRAY_BUFFER, b);\n          gl.bufferData(gl.ARRAY_BUFFER, new Float32Array(vtx), gl.STATIC_DRAW);\n          gl.vertexAttribPointer(0, 3, gl.FLOAT, false, 0, 0); gl.enableVertexAttribArray(0);\n          return {r,b}\n        },\n      // Destroy Vertex array and delete buffers.\n        destroyVA=function(va) {\n          if (va.b) gl.deleteBuffer(va.b); if (va.r) gl.va.deleteVertexArrayOES(va.r);\n        }\n      // Drawing function  \n        var M=[1,0,0,0,0,1,0,0,0,0,1,0,0,0,5,1];\n        var draw=function(p, tp, vtx, color, color2, ratio, texc, va, b,color3,r,g){\n            gl.useProgram(p); gl.uniformMatrix4fv(gl.getUniformLocation(p, \"mv\"),false,M);\n            gl.uniformMatrix4fv(gl.getUniformLocation(p, \"p\"),false, [5,0,0,0,0,5*(ratio||1),0,0,0,0,1,2,0,0,-1,0])\n            gl.uniform3fv(gl.getUniformLocation(p, \"color\"),new Float32Array(color));\n            gl.uniform3fv(gl.getUniformLocation(p, \"color2\"),new Float32Array(color2));\n            if (color3) gl.uniform3fv(gl.getUniformLocation(p, \"color3\"),new Float32Array(color3));\n            if (b) gl.uniform1fv(gl.getUniformLocation(p, \"b\"),(new Float32Array(counts[g])).map((x,i)=>b[g][i]||0));\n            if (texc) gl.uniform1i(gl.getUniformLocation(p, \"texc\"),0);\n            if (r) gl.uniform1f(gl.getUniformLocation(p,\"ratio\"),r);\n            var v; if (!va) v = createVA(vtx); else gl.va.bindVertexArrayOESOES(va.r);\n            gl.drawArrays(tp, 0, (va&&va.tcount)||vtx.length/3);\n            if (v) destroyVA(v);\n        }\n      // Compile the OPNS renderer. (sphere tracing)\n        var programs = [], genprog = grade=>compile(`${gl2?\"#version 300 es\":\"\"}\n             ${gl2?\"in\":\"attribute\"} vec4 position; ${gl2?\"out\":\"varying\"} vec4 Pos; uniform mat4 mv; uniform mat4 p;\n             void main() { Pos=mv*position; gl_Position = p*Pos; }`,\n            `${!gl2?\"#extension GL_EXT_frag_depth : enable\":\"#version 300 es\"}\n             precision highp float;  \n             uniform vec3 color; uniform vec3 color2; \n             uniform vec3 color3; uniform float b[${counts[grade]}];\n             uniform float ratio; ${gl2?\"out vec4 col;\":\"\"}\n             ${gl2?\"in\":\"varying\"} vec4 Pos; \n             float dist (in float z, in float y, in float x, in float[${counts[grade]}] b) {\n                ${this.nVector(1,[]).OPNS_GLSL(this.nVector(grade,[]), options.up)}\n                return ${grade!=tot-1?\"sign(sum)*sqrt(abs(sum))\":\"res\"};\n             }\n             vec3 trace_depth (in vec3 start, vec3 dir, in float thresh) {\n                vec3 orig=start; float lastd = 1000.0; const int count=${(options.maxSteps||64)};\n                float s =  sign(dist(start[0],start[1],start[2],b));\n                for (int i=0; i<count; i++) {\n                  float d = s*dist(start[0],start[1],start[2],b);\n                  if (d < thresh) return start - lastd*${(options.stepSize||0.25)}*dir*(thresh-d)/(lastd-d);\n                  lastd = d; start += dir*${(options.stepSize||0.25)}*d;\n                }\n                return orig;\n             }\n             void main() { \n               vec3 p = -5.0*normalize(color2); \n               vec3 dir = normalize((-Pos[0]/5.0)*color + color2 + vec3(0.0,Pos[1]/5.0*ratio,0.0));  p += 1.0*dir;\n               vec3 L = 5.0*normalize( -0.5*color + 0.85*color2 + vec3(0.0,-0.5,0.0) );\n               vec3 d2 = trace_depth( p , dir, ${grade!=tot-1?(options.thresh||0.2):\"0.0075\"} );\n               float dl2 = dot(d2-p,d2-p); const float h=0.1; \n               if (dl2>0.0) {\n                 vec3 n = normalize(vec3(\n                        dist(d2[0]+h,d2[1],d2[2],b)-dist(d2[0]-h,d2[1],d2[2],b),\n                        dist(d2[0],d2[1]+h,d2[2],b)-dist(d2[0],d2[1]-h,d2[2],b),\n                        dist(d2[0],d2[1],d2[2]+h,b)-dist(d2[0],d2[1],d2[2]-h,b)\n                      ));\n                 ${gl2?\"gl_FragDepth\":\"gl_FragDepthEXT\"} = dl2/50.0;\n                 ${gl2?\"col\":\"gl_FragColor\"} = vec4(max(0.2,abs(dot(n,normalize(L-d2))))*color3 + pow(abs(dot(n,normalize(normalize(L-d2)+dir))),100.0),1.0);\n               } else discard; \n             }`),genprog2D = grade=>compile(`${gl2?\"#version 300 es\":\"\"}\n             ${gl2?\"in\":\"attribute\"} vec4 position; ${gl2?\"out\":\"varying\"} vec4 Pos; uniform mat4 mv; uniform mat4 p;\n             void main() { Pos=mv*position; gl_Position = p*Pos; }`,\n            `${!gl2?\"#extension GL_EXT_frag_depth : enable\":\"#version 300 es\"}\n             precision highp float;  \n             uniform vec3 color; uniform vec3 color2; \n             uniform vec3 color3; uniform float b[${counts[grade]}];\n             uniform float ratio; ${gl2?\"out vec4 col;\":\"\"}\n             ${gl2?\"in\":\"varying\"} vec4 Pos; \n             float dist (in float z, in float y, in float x, in float[${counts[grade]}] b) {\n                ${this.nVector(1,[]).OPNS_GLSL(this.nVector(grade,[]), options.up)}\n                return ${grade!=tot-1?\"sqrt(abs(sum))\":\"res\"};\n             }\n             float trace_depth (in vec3 start, vec3 dir, in float thresh) {\n                vec3 orig=start; float lastd = 1000.0; const int count=${(options.maxSteps||64)};\n                float s = dist(start[0]*5.0,start[1]*5.0,start[2]*5.0,b);\n                s=s*s;\n                return 1.0-s*150.0;\n             }\n             void main() { \n               vec3 p = -5.0*normalize(color2); \n               vec3 dir = normalize((-Pos[0]/5.0)*color + color2 + vec3(0.0,Pos[1]/5.0*ratio,0.0));  p += 1.0*dir;\n               vec3 L = 5.0*normalize( -0.5*color + 0.85*color2 + vec3(0.0,-0.5,0.0) );\n               float d2 = trace_depth( p , dir, ${grade!=tot-1?(options.thresh||0.2):\"0.0075\"} );\n               if (d2>0.0) {\n                 ${gl2?\"gl_FragDepth\":\"gl_FragDepthEXT\"} = d2/50.0;\n                 ${gl2?\"col\":\"gl_FragColor\"} = vec4(d2*color3,d2);\n               } else discard;  \n             }`)\n      // canvas update will (re)render the content.            \n        var armed=0;\n        canvas.update = (x)=>{\n        // Start by updating canvas size if needed and viewport.\n          var s = getComputedStyle(canvas); if (s.width) { canvas.width = parseFloat(s.width); canvas.height = parseFloat(s.height); }\n          gl.viewport(0,0, canvas.width|0,canvas.height|0); var r=canvas.width/canvas.height;\n        // Defaults, resolve function input  \n          var a,p=[],l=[],t=[],c=[.5,.5,.5],alpha=0,lastpos=[-2,2,0.2]; gl.clear(gl.COLOR_BUFFER_BIT+gl.DEPTH_BUFFER_BIT); while (x.call) x=x();\n        // Loop over all items to render.  \n          for (var i=0,ll=x.length;i<ll;i++) { \n            var e=x[i]; while (e&&e.call) e=e(); if (e==undefined) continue;\n            if (typeof e == \"number\") { alpha=((e>>>24)&0xff)/255; c[0]=((e>>>16)&0xff)/255; c[1]=((e>>>8)&0xff)/255; c[2]=(e&0xff)/255; }\n            if (e instanceof Element){\n              var tt = options.spin?-performance.now()*options.spin/1000:-options.h||0; tt+=Math.PI/2; var r = canvas.height/canvas.width;\n              var g=tot-1; while(!e[g]&&g>1) g--;\n              if (!programs[tot-1-g]) programs[tot-1-g] = (options.up.find(x=>x.match&&x.match(\"z\")))?genprog(g):genprog2D(g);\n              gl.enable(gl.BLEND); gl.blendFunc(gl.ONE, gl.ONE_MINUS_SRC_ALPHA);\n              draw(programs[tot-1-g],gl.TRIANGLES,[-2,-2,0,-2,2,0,2,-2,0,-2,2,0,2,-2,0,2,2,0],[Math.cos(tt),0,-Math.sin(tt)],[Math.sin(tt),0,Math.cos(tt)],undefined,undefined,undefined,e,c,r,g);\n              gl.disable(gl.BLEND);\n            }\n          }\n          // if we're no longer in the page .. stop doing the work.\n          armed++; if (document.body.contains(canvas)) armed=0; if (armed==2) return;\n          canvas.value=x; if (options&&!options.animate) canvas.dispatchEvent(new CustomEvent('input'));\n          if (options&&options.animate) { requestAnimationFrame(canvas.update.bind(canvas,f,options)); }\n          if (options&&options.still) { canvas.value=x; canvas.dispatchEvent(new CustomEvent('input')); canvas.im.width=canvas.width; canvas.im.height=canvas.height; canvas.im.src = canvas.toDataURL(); }\n        }\n        // Basic mouse interactivity. needs more love.\n        var sel=-1; canvas.oncontextmenu = canvas.onmousedown = (e)=>{ e.preventDefault(); e.stopPropagation(); sel=-2;\n          var rc = canvas.getBoundingClientRect(), mx=(e.x-rc.left)/(rc.right-rc.left)*2-1, my=((e.y-rc.top)/(rc.bottom-rc.top)*-4+2)*canvas.height/canvas.width;\n          canvas.onwheel=e=>{e.preventDefault(); e.stopPropagation(); options.z = (options.z||5)+e.deltaY/100; if (!options.animate) requestAnimationFrame(canvas.update.bind(canvas,f,options));}\n          canvas.onmouseup=e=>sel=-1; canvas.onmouseleave=e=>sel=-1;\n          canvas.onmousemove=(e)=>{ \n            var rc = canvas.getBoundingClientRect(); \n            var mx =(e.movementX)/(rc.right-rc.left)*2, my=((e.movementY)/(rc.bottom-rc.top)*-2)*canvas.height/canvas.width;\n            if (sel==-2) { options.h =  (options.h||0)+mx; if (!options.animate) requestAnimationFrame(canvas.update.bind(canvas,f,options)); return; }; if (sel < 0) return;\n          }\n        }\n        canvas.value = f.call?f():f; canvas.options = options;\n        if (options&&options.still) {\n          var i=new Image(); canvas.im = i; return requestAnimationFrame(canvas.update.bind(canvas,f,options)),i;\n        } else return requestAnimationFrame(canvas.update.bind(canvas,f,options)),canvas;\n        \n      }\n    \n      \n    // webGL Graphing function. (for parametric defined objects)\n      static graphGL(f,options) {\n      // Create a canvas, webgl2 context and set some default GL options.\n        var canvas=document.createElement('canvas'); canvas.style.width=options.width||''; canvas.style.height=options.height||''; canvas.style.backgroundColor='#EEE';\n        if (options.width && options.width.match && options.width.match(/px/i)) canvas.width = parseFloat(options.width); if (options.height && options.height.match && options.height.match(/px/i)) canvas.height = parseFloat(options.height);\n        var gl=canvas.getContext('webgl',{alpha:options.alpha||false,antialias:true,preserveDrawingBuffer:options.still||true,powerPreference:'high-performance'}); \n        gl.enable(gl.DEPTH_TEST); gl.depthFunc(gl.LEQUAL); if (!options.alpha) gl.clearColor(240/255,240/255,240/255,1.0); gl.getExtension(\"OES_standard_derivatives\"); gl.va=gl.getExtension(\"OES_vertex_array_object\");\n      // Compile vertex and fragment shader, return program.  \n        var compile=(vs,fs)=>{ \n          var s=[gl.VERTEX_SHADER,gl.FRAGMENT_SHADER].map((t,i)=>{\n            var r=gl.createShader(t); gl.shaderSource(r,[vs,fs][i]); gl.compileShader(r);\n            return gl.getShaderParameter(r, gl.COMPILE_STATUS)&&r||console.error(gl.getShaderInfoLog(r));\n          });\n          var p = gl.createProgram(); gl.attachShader(p, s[0]); gl.attachShader(p, s[1]); gl.linkProgram(p);\n          gl.getProgramParameter(p, gl.LINK_STATUS)||console.error(gl.getProgramInfoLog(p));\n          return p;\n        };\n      // Create vertex array and buffers, upload vertices and optionally texture coordinates.  \n        var createVA=function(vtx, texc, idx, clr) {\n              var r = gl.va.createVertexArrayOES(); gl.va.bindVertexArrayOES(r);\n              var b = gl.createBuffer(); gl.bindBuffer(gl.ARRAY_BUFFER, b); \n              gl.bufferData(gl.ARRAY_BUFFER, new Float32Array(vtx), gl.STATIC_DRAW);\n              gl.vertexAttribPointer(0, 3, gl.FLOAT, false, 0, 0); gl.enableVertexAttribArray(0);\n              if (texc){\n                var b2=gl.createBuffer(); gl.bindBuffer(gl.ARRAY_BUFFER, b2);\n                gl.bufferData(gl.ARRAY_BUFFER, new Float32Array(texc), gl.STATIC_DRAW);\n                gl.vertexAttribPointer(1, 2, gl.FLOAT, false, 0, 0); gl.enableVertexAttribArray(1);\n              }\n              if (clr){\n                var b3=gl.createBuffer(); gl.bindBuffer(gl.ARRAY_BUFFER, b3);\n                gl.bufferData(gl.ARRAY_BUFFER, new Float32Array(clr), gl.STATIC_DRAW);\n                gl.vertexAttribPointer(1, 3, gl.FLOAT, false, 0, 0); gl.enableVertexAttribArray(1);\n              }\n              if (idx) {\n                var b4=gl.createBuffer(); gl.bindBuffer(gl.ELEMENT_ARRAY_BUFFER, b4);\n                gl.bufferData(gl.ELEMENT_ARRAY_BUFFER, new Uint16Array(idx), gl.STATIC_DRAW);\n              }\n              return {r,b,b2,b4,b3}\n            },\n      // Destroy Vertex array and delete buffers.\n            destroyVA=function(va) {\n              [va.b,va.b2,va.b4,va.b3].forEach(x=>{if(x) gl.deleteBuffer(x)}); if (va.r) gl.va.deleteVertexArrayOES(va.r);\n            }\n      // Default modelview matrix, convert camera to matrix (biquaternion->matrix)      \n        var M=[1,0,0,0,0,1,0,0,0,0,1,0,0,0,5,1], mtx = x=>{ var t=options.spin?performance.now()*options.spin/1000:options.h||0, t2=options.p||0;\n          var ct = Math.cos(t), st= Math.sin(t), ct2 = Math.cos(t2), st2 = Math.sin(t2), xx=options.posx||0, y=options.posy||0, z=options.posz||0, zoom=options.z||5;\n          if (tot==5) return [ct,st*-st2,st*ct2,0,0,ct2,st2,0,-st,ct*-st2,ct*ct2,0,xx*ct+z*-st,y*ct2+(xx*st+z*ct)*-st2,y*st2+xx*st+z*ct*ct2+zoom,1];\n          x=x.Normalized; var y=x.Mul(x.Dual),X=-x.e23,Y=-x.e13,Z=x.e12,W=x.s,m=Array(16);\n          var xx = X*X, xy = X*Y, xz = X*Z, xw = X*W, yy = Y*Y, yz = Y*Z, yw = Y*W, zz = Z*Z, zw = Z*W;\n          return [ 1-2*(yy+zz), 2*(xy+zw), 2*(xz-yw), 0, 2*(xy-zw), 1-2*(xx+zz), 2*(yz+xw), 0, 2*(xz+yw), 2*(yz-xw), 1-2*(xx+yy), 0, -2*y.e23, -2*y.e13, 2*y.e12+5, 1];\n        }\n      // Render the given vertices. (autocreates/destroys vertex array if not supplied).  \n        var draw=function(p, tp, vtx, color, color2, ratio, texc, va){\n          gl.useProgram(p); gl.uniformMatrix4fv(gl.getUniformLocation(p, \"mv\"),false,M); \n          gl.uniformMatrix4fv(gl.getUniformLocation(p, \"p\"),false, [5,0,0,0,0,5*(ratio||2),0,0,0,0,1,2,0,0,-1,0])\n          gl.uniform3fv(gl.getUniformLocation(p, \"color\"),new Float32Array(color));\n          gl.uniform3fv(gl.getUniformLocation(p, \"color2\"),new Float32Array(color2));\n          if (texc) gl.uniform1i(gl.getUniformLocation(p, \"texc\"),0);\n          var v; if (!va) v = createVA(vtx, texc); else gl.va.bindVertexArrayOES(va.r);\n          if (va && va.b4) {\n            gl.drawElements(tp, va.tcount, gl.UNSIGNED_SHORT, 0);\n          } else {\n            gl.drawArrays(tp, 0, (va&&va.tcount)||vtx.length/3);\n          }  \n          if (v) destroyVA(v);\n        }\n      // Program for the geometry. Derivative based normals. Basic lambert shading.    \n        var program = compile(`attribute vec4 position; varying vec4 Pos; uniform mat4 mv; uniform mat4 p; \n                 void main() { gl_PointSize=6.0; Pos=mv*position; gl_Position = p*Pos; }`,\n                `#extension GL_OES_standard_derivatives : enable\n                 precision highp float; uniform vec3 color; uniform vec3 color2; varying vec4 Pos; \n                 void main() { vec3 ldir = normalize(Pos.xyz - vec3(1.0,1.0,2.0));\n                 vec3 normal = normalize(cross(dFdx(Pos.xyz), dFdy(Pos.xyz))); float l=dot(normal,ldir);\n                 vec3 E = normalize(-Pos.xyz); vec3 R = normalize(reflect(ldir,normal));  \n                 gl_FragColor = vec4(max(0.0,l)*color+vec3(0.5*pow(max(dot(R,E),0.0),20.0))+color2, 1.0);  }`);\n        var programcol = compile(`attribute vec4 position; attribute vec3 col; varying vec3 Col; varying vec4 Pos; uniform mat4 mv; uniform mat4 p; \n                 void main() { gl_PointSize=6.0; Pos=mv*position; gl_Position = p*Pos; Col=col; }`,\n                `#extension GL_OES_standard_derivatives : enable\n                 precision highp float; uniform vec3 color; uniform vec3 color2; varying vec4 Pos; varying vec3 Col; \n                 void main() { vec3 ldir = normalize(Pos.xyz - vec3(1.0,1.0,2.0));\n                 vec3 normal = normalize(cross(dFdx(Pos.xyz), dFdy(Pos.xyz))); float l=dot(normal,ldir);\n                 vec3 E = normalize(-Pos.xyz); vec3 R = normalize(reflect(ldir,normal));  \n                 gl_FragColor = vec4(max(0.3,l)*Col+vec3(0.5*pow(max(dot(R,E),0.0),20.0))+color2, 1.0);  }`);\n      // Create a font texture, lucida console or otherwise monospaced.\n        var fw=22, font = Object.assign(document.createElement('canvas'),{width:94*fw,height:32}), \n            ctx = Object.assign(font.getContext('2d'),{font:'bold 32px lucida console, monospace'}),\n            ftx = gl.createTexture(); gl.activeTexture(gl.TEXTURE0); gl.bindTexture(gl.TEXTURE_2D, ftx);\n            for (var i=33; i<127; i++) ctx.fillText(String.fromCharCode(i),(i-33)*fw,26);\n            // 2.0 gl.texImage2D(gl.TEXTURE_2D,0,gl.RGBA,94*fw,32,0,gl.RGBA,gl.UNSIGNED_BYTE,font);\n            gl.texImage2D(gl.TEXTURE_2D, 0, gl.RGBA, gl.RGBA, gl.UNSIGNED_BYTE, font);\n\n            gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.LINEAR); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.LINEAR);\n            gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.CLAMP_TO_EDGE); gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_T, gl.CLAMP_TO_EDGE);  \n      // Font rendering program. Renders billboarded fonts, transforms offset passed as color2.\n        var program2 = compile(`attribute vec4 position; attribute vec2 texc; varying vec2 tex; varying vec4 Pos; uniform mat4 mv; uniform mat4 p; uniform vec3 color2; \n                 void main() { tex=texc; gl_PointSize=6.0; vec4 o=mv*vec4(color2,0.0); Pos=(-1.0/(o.z-mv[3][2]))*position+vec4(mv[3][0],mv[3][1],mv[3][2],0.0)+o; gl_Position = p*Pos; }`,\n                `precision highp float; uniform vec3 color; varying vec4 Pos; varying vec2 tex; \n                 uniform sampler2D texm; void main() { vec4 c = texture2D(texm,tex); if (c.a<0.01) discard; gl_FragColor = vec4(color,c.a);}`);\n      // Conformal space needs a bit extra magic to extract euclidean parametric representations.\n        if (tot==5 && options.conformal) var ninf = Element.Coeff(4,1).Add(Element.Coeff(5,1)), no = Element.Coeff(4,0.5).Sub(Element.Coeff(5,0.5));\n        var interprete = (x)=>{\n          if (!(x instanceof Element)) return { tp:0 };\n          // tp = { 0:unknown 1:point 2:line, 3:plane, 4:circle, 5:sphere\n          var X2 = (x.Mul(x)).s, tp=0, weight2, opnix = ninf.Wedge(x), ipnix = ninf.LDot(x), \n              attitude, pos, normal, tg,btg,epsilon = 0.001/(options.scale||1), I3=Element.Coeff(16,-1);\n          var x2zero = Math.abs(X2) < epsilon, ipnixzero = ipnix.VLength < epsilon, opnixzero = opnix.VLength < epsilon;\n          if (opnixzero && ipnixzero) {                 // free flat\n          } else if (opnixzero && !ipnixzero) {         // bound flat (lines)\n            attitude = no.Wedge(ninf).LDot(x); \n            weight2 = Math.abs(attitude.LDot(attitude).s)**.5;\n            pos = attitude.LDot(x.Reverse); //Inverse);\n            pos = [-pos.e15/pos.e45,-pos.e25/pos.e45,-pos.e34/pos.e45];\n            if (x.Grade(3).VLength) {\n              normal = [attitude.e1/weight2,attitude.e2/weight2,attitude.e3/weight2]; tp=2; \n            } else {\n              normal = Element.LDot(Element.Mul(attitude,1/weight2),I3).Normalized;\n              var r=normal.Mul(Element.Coeff(3,1)); if (r[0]==-1) r[0]=1; else {r[0]+=1; r=r.Normalized;}\n              tg = [...r.Mul(Element.Coeff(1,1)).Mul(r.Conjugate)].slice(1,4);\n              btg = [...r.Mul(Element.Coeff(2,1)).Mul(r.Conjugate)].slice(1,4);\n              normal = [...normal.slice(1,4)]; tp=3;\n            }\n          } else if (!opnixzero && ipnixzero) {         // dual bound flat\n          } else if (x2zero) {                          // bound vec,biv,tri (points)\n            attitude = ninf.Wedge(no).LDot(ninf.Wedge(x)); \n            pos = [...(Element.LDot(1/(ninf.LDot(x)).s,x)).slice(1,4)].map(x=>-x);\n            tp=1; \n          } else if (!x2zero) {                          // round (point pair,circle,sphere)\n            tp = x.Grade(3).VLength?4:x.Grade(2).VLength?6:5; \n            var nix  = ninf.Wedge(x), nix2 = (nix.Mul(nix)).s;\n            attitude = ninf.Wedge(no).LDot(nix);\n            pos = [...(x.Mul(ninf).Mul(x)).slice(1,4)].map(x=>-x/(2.0*nix2));\n            weight2 = Math.abs((x.LDot(x)).s / nix2)**.5;\n            if (tp==4) {\n              if (x.LDot(x).s < 0) { weight2 = -weight2; }\n              normal = Element.LDot(Element.Mul(attitude,1/weight2),I3).Normalized;\n              var r=normal.Mul(Element.Coeff(3,1)); if (r[0]==-1) r[0]=1; else {r[0]+=1; r=r.Normalized;}\n              tg = [...r.Mul(Element.Coeff(1,1)).Mul(r.Conjugate)].slice(1,4);\n              btg = [...r.Mul(Element.Coeff(2,1)).Mul(r.Conjugate)].slice(1,4);\n              normal = [...normal.slice(1,4)]; \n            } else if (tp==6) {\n              weight2 = (x.LDot(x).s < 0)?-(weight2):weight2;\n              normal = Element.Mul(attitude.Normalized,weight2).slice(1,4);\n            } else {\n              normal = [...((Element.LDot(Element.Mul(attitude,1/weight2),I3)).Normalized).slice(1,4)];\n            }\n          }\n          return {tp,pos:pos?pos.map(x=>x*(options.scale||1)):[0,0,0],normal,tg,btg,weight2:weight2*(options.scale||1)}\n        };                 \n      // canvas update will (re)render the content.            \n        var armed=0,sphere,e14 = Element.Coeff(14,1);\n        canvas.update = (x)=>{\n        // restore from still..\n          if (options && !options.still && canvas.im && canvas.im.parentElement) { canvas.im.parentElement.insertBefore(canvas,canvas.im); canvas.im.parentElement.removeChild(canvas.im); }\n        // Start by updating canvas size if needed and viewport.\n          var s = getComputedStyle(canvas); if (s.width) { canvas.width = parseFloat(s.width); canvas.height = parseFloat(s.height); }\n          gl.viewport(0,0, canvas.width|0,canvas.height|0); var r=canvas.width/canvas.height;\n        // Defaults, resolve function input  \n          var a,p=[],l=[],t=[],c=[.5,.5,.5],alpha=0,lastpos=[-2,2,0.2]; gl.clear(gl.COLOR_BUFFER_BIT+gl.DEPTH_BUFFER_BIT); while (x.call) x=x();\n        // Create default camera matrix and initial lastposition (contra-compensated for camera)  \n          M = mtx(options.camera); lastpos = options.camera.Normalized.Conjugate.Mul(((a=new this()).set(lastpos,11),a)).Mul(options.camera.Normalized).slice(11,14);\n        // Grid.\n          if (options.grid) {\n            if (!options.gridLines) { options.gridLines=[[],[],[]]; for (var i=-5; i<=5; i++) {\n                options.gridLines[0].push(i,0,5, i,0,-5, 5,0,i, -5,0,i); options.gridLines[1].push(i,5,0, i,-5,0, 5,i,0, -5,i,0); options.gridLines[2].push(0,i,5, 0,i,-5, 0,5,i, 0,-5,i);\n            }}\n            gl.depthMask(false);\n              draw(program,gl.LINES,options.gridLines[0],[0,0,0],[.6,1,.6],r); draw(program,gl.LINES,options.gridLines[1],[0,0,0],[1,.8,.8],r); draw(program,gl.LINES,options.gridLines[2],[0,0,0],[.8,.8,1],r);\n            gl.depthMask(true);\n          }\n        // Z-buffer override.\n          if (options.noZ) gl.depthMask(false);  \n        // Loop over all items to render.  \n          for (var i=0,ll=x.length;i<ll;i++) { \n            var e=x[i]; while (e&&e.call&&e.length==0) e=e(); if (e==undefined) continue;\n          // CGA\n            if (tot==5 && options.conformal) {\n              if (e instanceof Array && e.length==2) { e.forEach(x=>{ while (x.call) x=x.call(); x=interprete(x);l.push.apply(l,x.pos); });  var d = {tp:-1}; }\n              else if (e instanceof Array && e.length==3) { e.forEach(x=>{ while (x.call) x=x.call(); x=interprete(x);t.push.apply(t,x.pos); });  var d = {tp:-1}; }\n              else var d = interprete(e);\n              if (d.tp) lastpos=d.pos;\n              if (d.tp==1) p.push.apply(p,d.pos);\n              if (d.tp==2) { l.push.apply(l,d.pos.map((x,i)=>x-d.normal[i]*10)); l.push.apply(l,d.pos.map((x,i)=>x+d.normal[i]*10)); }\n              if (d.tp==3) { t.push.apply(t,d.pos.map((x,i)=>x+d.tg[i]+d.btg[i])); t.push.apply(t,d.pos.map((x,i)=>x-d.tg[i]+d.btg[i])); t.push.apply(t,d.pos.map((x,i)=>x+d.tg[i]-d.btg[i])); \n                             t.push.apply(t,d.pos.map((x,i)=>x-d.tg[i]+d.btg[i])); t.push.apply(t,d.pos.map((x,i)=>x+d.tg[i]-d.btg[i])); t.push.apply(t,d.pos.map((x,i)=>x-d.tg[i]-d.btg[i])); }\n              if (d.tp==4) {\n                var ne=0,la=0;\n                if (d.weight2<0) { c[0]=1;c[1]=0;c[2]=0; }\n                for (var j=0; j<65; j++) {\n                  ne = d.pos.map((x,i)=>x+Math.cos(j/32*Math.PI)*d.weight2*d.tg[i]+Math.sin(j/32*Math.PI)*d.weight2*d.btg[i]); if (ne&&la&&(d.weight2>0||j%2==0)) { l.push.apply(l,la); l.push.apply(l,ne); }; la=ne;\n                }\n              }               \n              if (d.tp==6) {\n                if (d.weight2<0) { c[0]=1;c[1]=0;c[2]=0; }\n                if (options.useUnnaturalLineDisplayForPointPairs) {\n                  l.push.apply(l,d.pos.map((x,i)=>x-d.normal[i]*(options.scale||1)));\n                  l.push.apply(l,d.pos.map((x,i)=>x+d.normal[i]*(options.scale||1)));\n                } \n                p.push.apply(p,d.pos.map((x,i)=>x-d.normal[i]*(options.scale||1)));\n                p.push.apply(p,d.pos.map((x,i)=>x+d.normal[i]*(options.scale||1)));\n              }\n              if (d.tp==5) {\n                if (!sphere) {\n                  var pnts = [], tris=[], S=Math.sin, C=Math.cos, pi=Math.PI, W=96, H=48;\n                  for (var j=0; j<W+1; j++) for (var k=0; k<H; k++) {\n                    pnts.push( [S(2*pi*j/W)*S(pi*k/(H-1)), C(2*pi*j/W)*S(pi*k/(H-1)), C(pi*k/(H-1))]);\n                    if (j && k) {\n                      tris.push.apply(tris, pnts[(j-1)*H+k-1]);tris.push.apply(tris, pnts[(j-1)*H+k]);tris.push.apply(tris, pnts[j*H+k-1]);\n                      tris.push.apply(tris, pnts[j*H+k-1]); tris.push.apply(tris, pnts[(j-1)*H+k]); tris.push.apply(tris, pnts[j*H+k]);\n                  }}\n                  sphere = { va : createVA(tris,undefined) }; sphere.va.tcount = tris.length/3;\n                }\n                var oldM = M;\n                M=[].concat.apply([],Element.Mul([[d.weight2,0,0,0],[0,d.weight2,0,0],[0,0,d.weight2,0],[d.pos[0],d.pos[1],d.pos[2],1]],[[M[0],M[1],M[2],M[3]],[M[4],M[5],M[6],M[7]],[M[8],M[9],M[10],M[11]],[M[12],M[13],M[14],M[15]]])).map(x=>x.s);\n                gl.enable(gl.BLEND); gl.blendFunc(gl.CONSTANT_ALPHA, gl.ONE_MINUS_CONSTANT_ALPHA); gl.blendColor(1,1,1,0.5); gl.enable(gl.CULL_FACE)\n                draw(program,gl.TRIANGLES,undefined,c,[0,0,0],r,undefined,sphere.va);\n                gl.disable(gl.BLEND); gl.disable(gl.CULL_FACE);\n                M = oldM;\n              }\n              if (i==ll-1 || d.tp==0) {\n              // render triangles, lines, points.\n                if (alpha) { gl.enable(gl.BLEND); gl.blendFunc(gl.CONSTANT_ALPHA, gl.ONE_MINUS_CONSTANT_ALPHA); gl.blendColor(1,1,1,1-alpha); }\n                if (t.length) { draw(program,gl.TRIANGLES,t,c,[0,0,0],r); t.forEach((x,i)=>{ if (i%9==0) lastpos=[0,0,0]; lastpos[i%3]+=x/3; }); t=[];  }\n                if (l.length) { draw(program,gl.LINES,l,[0,0,0],c,r); var l2=l.length-1; lastpos=[(l[l2-2]+l[l2-5])/2,(l[l2-1]+l[l2-4])/2+0.1,(l[l2]+l[l2-3])/2]; l=[]; }\n                if (p.length) { draw(program,gl.POINTS,p,[0,0,0],c,r); lastpos = p.slice(-3); lastpos[0]-=0.075; lastpos[1]+=0.075; p=[]; }\n                if (alpha) gl.disable(gl.BLEND);\n                // we could also be an object with cached vertex array of triangles ..   \n                  if (e instanceof Object && e.data) {\n                    // Create the vertex array and store it for re-use.\n                    if (!e.va) {\n                      var et=[],et2=[],et3=[],lc=0,pc=0,tc=0; e.data.forEach(e=>{ \n                        if (e instanceof Array && e.length==3) { tc++; e.forEach(x=>{ while (x.call) x=x.call(); x=interprete(x);et3.push.apply(et3,x.pos); });  var d = {tp:-1}; }\n                        else {\n                          var d = interprete(e); \n                          if (d.tp==1) { pc++; et.push(...d.pos);  }\n                          if (d.tp==2) { lc++; et2.push(...d.pos.map((x,i)=>x-d.normal[i]*10),...d.pos.map((x,i)=>x+d.normal[i]*10)); }\n                        }\n                      });\n                      e.va = createVA(et,undefined); e.va.tcount = pc;\n                      e.va2 = createVA(et2,undefined); e.va2.tcount = lc*2;\n                      e.va3 = createVA(et3,undefined); e.va3.tcount = tc*3;\n                    }\n                    // render the vertex array.\n                    if (e.va.tcount) draw(program,gl.POINTS,undefined,[0,0,0],c,r,undefined,e.va);\n                    if (e.va2.tcount) draw(program,gl.LINES,undefined,[0,0,0],c,r,undefined,e.va2);\n                    if (e.va3.tcount) draw(program,gl.TRIANGLES,undefined,[0,0,0],c,r,undefined,e.va3);\n                  }\n              // setup a new color  \n                if (typeof e == \"number\") { alpha=((e>>>24)&0xff)/255; c[0]=((e>>>16)&0xff)/255; c[1]=((e>>>8)&0xff)/255; c[2]=(e&0xff)/255; }\n                if (typeof(e)=='string') {\n                  gl.enable(gl.BLEND); gl.blendFunc(gl.SRC_ALPHA,gl.ONE_MINUS_SRC_ALPHA); \n                  draw(program2,gl.TRIANGLES, \n                       [...Array(e.length*6*3)].map((x,i)=>{ var x=0,z=-0.2, o=x+(i/18|0)*1.1; return (0.05*(options.z||5))*[o,-1,z,o+1.2,-1,z,o,1,z,o+1.2,-1,z,o+1.2,1,z,o,1,z][i%18]}),c,lastpos,r,\n                       [...Array(e.length*6*2)].map((x,i)=>{ var o=(e.charCodeAt(i/12|0)-33)/94; return [o,1,o+1/94,1,o,0,o+1/94,1,o+1/94,0,o,0][i%12]})); gl.disable(gl.BLEND); lastpos[1]-=0.18;\n                }\n              }\n              continue;\n            }\n          // PGA   \n          // Convert planes to polygons.\n            if (e instanceof Element && e.Grade(1).Length) {\n              var m = Element.Add(1, Element.Mul(e.Normalized, Element.Coeff(3,1))).Normalized, e0 = 0;\n              e=Element.sw(m,[[-1,-1],[-1,1],[1,1],[-1,-1],[1,1],[1,-1]].map(([x,z])=>Element.Trivector(x,e0,z,1)));\n            }\n          // Convert lines to line segments.  \n            if (e instanceof Element && e.Grade(2).Length) \n               e=[e.LDot(e14).Wedge(e).Add(e.Wedge(Element.Coeff(1,1)).Mul(Element.Coeff(0,-500))),e.LDot(e14).Wedge(e).Add(e.Wedge(Element.Coeff(1,1)).Mul(Element.Coeff(0,500)))];\n          // If euclidean point, store as point, store line segments and triangles.\n            if (e.e123) p.push.apply(p,e.slice(11,14).map((y,i)=>(i==0?1:-1)*y/e[14]).reverse());\n            if (e instanceof Array && e.length==2) l=l.concat.apply(l,e.map(x=>[...x.slice(11,14).map((y,i)=>(i==0?1:-1)*y/x[14]).reverse()])); \n            if (e instanceof Array && e.length%3==0) t=t.concat.apply(t,e.map(x=>[...x.slice(11,14).map((y,i)=>(i==0?1:-1)*y/x[14]).reverse()]));\n          // Render orbits of parametrised motors\n            function sw_mot_orig(A,R){\n              var a0=A[0],a1=A[5],a2=A[6],a3=A[7],a4=A[8],a5=A[9],a6=A[10],a7=A[15];\n              R[2] = -2*(a0*a3+a4*a7-a6*a2-a5*a1);\n              R[1] = -2*(a4*a1-a0*a2-a6*a3+a5*a7);\n              R[0] =  2*(a0*a1+a4*a2+a5*a3+a6*a7);\n              return R\n            }\n            if ( e.call && e.length==1) { var count=e.dx||64;\n              for (var xx,o=new Float32Array(3),ii=0; ii<count; ii++) {\n                if (ii>1) l.push(xx[0],xx[1],xx[2]);\n                xx = sw_mot_orig(e(ii/(count-1)),o); //Element.sw(e(ii/(count-1)),o);\n                l.push(xx[0],xx[1],xx[2]);\n              }\n            }\n            if ( e.call && e.length==2 && !e.va) { var countx=e.dx||64,county=e.dy||32;\n              var temp=new Float32Array(3*countx*county),o=new Float32Array(3),et=[];\n              for (var pp=0,ii=0; ii<countx; ii++) for (var jj=0; jj<county; jj++,pp+=3) temp.set(sw_mot_orig(e(ii/(countx-1),jj/(county-1)),o),pp);\n              for (ii=0; ii<countx-1; ii++) for (var jj=0; jj<county; jj++) et.push((ii+0)*county+(jj+0),(ii+0)*county+(jj+1),(ii+1)*county+(jj+1),(ii+0)*county+(jj+0),(ii+1)*county+(jj+1),(ii+1)*county+(jj+0));\n              e.va = createVA(temp,undefined,et.map(x=>x%(countx*county))); e.va.tcount = (countx-1)*county*2*3;\n            }\n          // we could also be an object with cached vertex array of triangles ..   \n            if (e.va || (e instanceof Object && e.data)) {\n              // Create the vertex array and store it for re-use.\n              if (!e.va) {\n                if (e.idx) {\n                  var et = e.data.map(x=>[...x.slice(11,14).map((y,i)=>(i==0?1:-1)*y/x[14]).reverse()]).flat();\n                } else {\n                  var et=[]; e.data.forEach(e=>{if (e instanceof Array && e.length==3) et=et.concat.apply(et,e.map(x=>[...x.slice(11,14).map((y,i)=>(i==0?1:-1)*y/x[14]).reverse()]));});\n                }\n                e.va = createVA(et,undefined,e.idx,e.color?new Float32Array(e.color):undefined); e.va.tcount = (e.idx && e.idx.length)?e.idx.length:e.data.length*3;\n              }\n              // render the vertex array.\n              if (e.transform) { M=mtx(options.camera.Mul(e.transform)); }\n              if (alpha) { gl.enable(gl.BLEND); gl.blendFunc(gl.CONSTANT_ALPHA, gl.ONE_MINUS_CONSTANT_ALPHA); gl.blendColor(1,1,1,1-alpha); }\n              draw(e.color?programcol:program,gl.TRIANGLES,t,c,[0,0,0],r,undefined,e.va);\n              if (alpha) gl.disable(gl.BLEND);\n              if (e.transform) { M=mtx(options.camera); }\n            }\n          // if we're a number (color), label or the last item, we output the collected items.  \n            else if (!isNaN(e) || i==ll-1 || typeof e == 'string') {\n            // render triangles, lines, points.\n              if (alpha) { gl.enable(gl.BLEND); gl.blendFunc(gl.CONSTANT_ALPHA, gl.ONE_MINUS_CONSTANT_ALPHA); gl.blendColor(1,1,1,1-alpha); }\n              if (t.length) { draw(program,gl.TRIANGLES,t,c,[0,0,0],r); t.forEach((x,i)=>{ if (i%9==0) lastpos=[0,0,0]; lastpos[i%3]+=x/3; }); t=[];  }\n              if (l.length) { draw(program,gl.LINES,l,[0,0,0],c,r); var l2=l.length-1; lastpos=[(l[l2-2]+l[l2-5])/2,(l[l2-1]+l[l2-4])/2+0.1,(l[l2]+l[l2-3])/2]; l=[]; }\n              if (p.length) { draw(program,gl.POINTS,p,[0,0,0],c,r); lastpos = p.slice(-3); lastpos[0]-=0.075; lastpos[1]+=0.075; p=[]; }\n              if (alpha) gl.disable(gl.BLEND);\n            // setup a new color  \n              if (typeof e == \"number\") { alpha=((e>>>24)&0xff)/255; c[0]=((e>>>16)&0xff)/255; c[1]=((e>>>8)&0xff)/255; c[2]=(e&0xff)/255; }\n            // render a label  \n              if (typeof(e)=='string') {\n                gl.enable(gl.BLEND); gl.blendFunc(gl.SRC_ALPHA,gl.ONE_MINUS_SRC_ALPHA); \n                draw(program2,gl.TRIANGLES, \n                     [...Array(e.length*6*3)].map((x,i)=>{ var x=0,z=-0.2, o=x+(i/18|0)*1.1; return 0.25*[o,-1,z,o+1.2,-1,z,o,1,z,o+1.2,-1,z,o+1.2,1,z,o,1,z][i%18]}),c,lastpos,r,\n                     [...Array(e.length*6*2)].map((x,i)=>{ var o=(e.charCodeAt(i/12|0)-33)/94; return [o,1,o+1/94,1,o,0,o+1/94,1,o+1/94,0,o,0][i%12]})); gl.disable(gl.BLEND); lastpos[1]-=0.18;\n              }\n            }  \n          }; \n          // if we're no longer in the page .. stop doing the work.\n          armed++; if (document.body.contains(canvas)) armed=0; if (armed==2) return;\n          canvas.value=x; if (options&&!options.animate) canvas.dispatchEvent(new CustomEvent('input')); canvas.options=options;\n          if (options&&options.animate) { requestAnimationFrame(canvas.update.bind(canvas,f,options)); }\n          if (options&&options.still) { canvas.value=x; canvas.dispatchEvent(new CustomEvent('input')); canvas.im.style.width=canvas.style.width; canvas.im.style.height=canvas.style.height; canvas.im.src = canvas.toDataURL(); \n            var p=canvas.parentElement;  if (p) { p.insertBefore(canvas.im,canvas); p.removeChild(canvas); }\n          }\n        }\n        // Basic mouse interactivity. needs more love.\n        var sel=-1; canvas.oncontextmenu = canvas.onmousedown = (e)=>{e.preventDefault(); e.stopPropagation();  if (e.detail===0) return; \n          var rc = canvas.getBoundingClientRect(), mx=(e.x-rc.left)/(rc.right-rc.left)*2-1, my=((e.y-rc.top)/(rc.bottom-rc.top)*-4+2)*canvas.height/canvas.width;\n          sel = (e.button==2)?-3:-2; canvas.value.forEach((x,i)=>{\n            x = interprete(x); if (x.tp==1) {\n              var pos2 = Element.Mul( [[M[0],M[4],M[8],M[12]],[M[1],M[5],M[9],M[13]],[M[2],M[6],M[10],M[14]],[M[3],M[7],M[11],M[15]]], [...x.pos,1]).map(x=>x.s);\n              pos2 = Element.Mul( [[5,0,0,0],[0,5*(r||2),0,0],[0,0,1,-1],[0,0,2,0]], pos2).map(x=>x.s).map((x,i,a)=>x/a[3]);\n              if ((mx-pos2[0])**2 + (my-pos2[1])**2 < 0.01) sel=i;\n            }\n          });\n          canvas.onwheel=e=>{e.preventDefault(); e.stopPropagation(); options.z = (options.z||5)+e.deltaY/100; if (!options.animate) requestAnimationFrame(canvas.update.bind(canvas,f,options));}\n          canvas.onmouseup=e=>sel=-1; canvas.onmouseleave=e=>sel=-1;\n          var tx,ty; canvas.ontouchstart = (e)=>{e.preventDefault();  canvas.focus(); var x = e.changedTouches[0].pageX, y = e.changedTouches[0].pageY; tx=x; ty=y; }\n          canvas.ontouchmove = function (e) { e.preventDefault();\n             var x = e.changedTouches[0].pageX, y = e.changedTouches[0].pageY, mx = (x-(tx||x))/1000, my = -(y-(ty||y))/1000; tx=x; ty=y; \n             options.h = (options.h||0)+mx; options.p = Math.max(-Math.PI/2,Math.min(Math.PI/2, (options.p||0)+my)); if (!options.animate) requestAnimationFrame(canvas.update.bind(canvas,f,options)); return;  \n          };\n          canvas.onmousemove=(e)=>{ \n            var rc = canvas.getBoundingClientRect(), x=interprete(canvas.value[sel]);\n            var mx =(e.movementX)/(rc.right-rc.left)*2, my=((e.movementY)/(rc.bottom-rc.top)*-2)*canvas.height/canvas.width;\n            if (sel==-2) { options.h =  (options.h||0)+mx; options.p = Math.max(-Math.PI/2,Math.min(Math.PI/2, (options.p||0)+my)); if (!options.animate) requestAnimationFrame(canvas.update.bind(canvas,f,options)); return; }; \n            if (sel==-3) { var ct = Math.cos(options.h||0), st= Math.sin(options.h||0), ct2 = Math.cos(options.p||0), st2 = Math.sin(options.p||0);\n              if (e.shiftKey) { options.posy = (options.posy||0)+my; } else { options.posx = (options.posx||0)+mx*ct+my*st; options.posz = (options.posz||0)+mx*-st+my*ct*ct2; } if (!options.animate) requestAnimationFrame(canvas.update.bind(canvas,f,options));return; }; if (sel < 0) return;\n            x.pos[0] += (e.buttons!=2)?Math.cos(-(options.h||0))*mx:Math.sin((options.h||0))*my; x.pos[1]+=(e.buttons!=2)?my:0; x.pos[2]+=(e.buttons!=2)?Math.sin(-(options.h||0))*mx:Math.cos((options.h||0))*my;\n            canvas.value[sel].set(Element.Mul(ninf,(x.pos[0]**2+x.pos[1]**2+x.pos[2]**2)*0.5).Sub(no)); canvas.value[sel].set(x.pos,1);\n            if (!options.animate) requestAnimationFrame(canvas.update.bind(canvas,f,options));\n          }\n        }\n        canvas.value = f.call?f():f; canvas.options=options;\n        if (options&&options.still) {\n          var i=new Image(); canvas.im = i; return requestAnimationFrame(canvas.update.bind(canvas,f,options)),canvas;\n        } else return requestAnimationFrame(canvas.update.bind(canvas,f,options)),canvas;\n      }  \n    \n    // The inline function is a js to js translator that adds operator overloading and algebraic literals.\n    // It can be called with a function, a string, or used as a template function.  \n      static inline(intxt) {\n      // If we are called as a template function. \n        if (arguments.length>1 || intxt instanceof Array) {\n          var args=[].slice.call(arguments,1);\n          return res.inline(new Function(args.map((x,i)=>'_template_'+i).join(),'return ('+intxt.map((x,i)=>(x||'')+(args[i]&&('_template_'+i)||'')).join('')+')')).apply(res,args);\n        }\n      // Get the source input text.    \n        var txt = (intxt instanceof Function)?intxt.toString():`function(){return (${intxt})}`;\n      // Our tokenizer reads the text token by token and stores it in the tok array (as type/token tuples).   \n        var tok = [], resi=[], t, tokens = [/^[\\s\\uFFFF]|^[\\u000A\\u000D\\u2028\\u2029]|^\\/\\/[^\\n]*\\n|^\\/\\*[\\s\\S]*?\\*\\//g,                 // 0: whitespace/comments\n          /^\\\"\\\"|^\\'\\'|^\\\".*?[^\\\\]\\\"|^\\'.*?[^\\\\]\\'|^\\`[\\s\\S]*?[^\\\\]\\`/g,                                                                // 1: literal strings\n          /^\\d+[.]{0,1}\\d*[ei][\\+\\-_]{0,1}\\d*|^\\.\\d+[ei][\\+\\-_]{0,1}\\d*|^e_\\d*/g,                                                       // 2: literal numbers in scientific notation (with small hack for i and e_ asciimath)\n          /^\\d+[.]{0,1}\\d*[E][+-]{0,1}\\d*|^\\.\\d+[E][+-]{0,1}\\d*|^0x\\d+|^\\d+[.]{0,1}\\d*|^\\.\\d+|^\\(\\/.*[^\\\\]\\/\\)/g,                       // 3: literal hex, nonsci numbers and regex (surround regex with extra brackets!)\n          /^(\\.Normalized|\\.Length|\\.\\.\\.|>>>=|===|!==|>>>|<<=|>>=|=>|\\|\\||[<>\\+\\-\\*%&|^\\/!\\=]=|\\*\\*|\\+\\+|\\-\\-|<<|>>|\\&\\&|\\^\\^|^[{}()\\[\\];.,<>\\+\\-\\*%|&^!~?:=\\/]{1})/g,   // 4: punctuator\n          /^[A-Za-z0-9_]*/g]                                                                                                            // 5: identifier\n        while (txt.length) for(t in tokens) if(resi=txt.match(tokens[t])){ tok.push([t|0,resi[0]]); txt=txt.slice(resi[0].length); break;} // tokenise \n      // Translate algebraic literals. (scientific e-notation to \"this.Coeff\"\n        tok=tok.map(t=>(t[0]==2)?[2,'Element.Coeff('+basis.indexOf('e'+(t[1].split(/e_|e|i/)[1]||1))+','+parseFloat(t[1][0]=='e'?1:t[1].split(/e_|e|i/)[0])+')']:t);\n      // We support two syntaxes, standard js or if you pass in a text, asciimath.       \n        var syntax = (intxt instanceof Function)?[[['.Normalized','Normalize',2],['.Length','Length',2]],[['~','Conjugate',1],['!','Dual',1]],[['**','Pow',0,1]],[['^','Wedge'],['&','Vee'],['<<','LDot']],[['*','Mul'],['/','Div']],[['|','Dot']],[['>>>','sw',0,1]],[['-','Sub'],['+','Add']],[['==','eq'],['!=','neq'],['<','lt'],['>','gt'],['<=','lte'],['>=','gte']]]\n                                                :[[['pi','Math.PI'],['sin','Math.sin']],[['ddot','this.Reverse'],['tilde','this.Involute'],['hat','this.Conjugate'],['bar','this.Dual']],[['^','Pow',0,1]],[['^^','Wedge'],['*','LDot']],[['**','Mul'],['/','Div']],[['-','Sub'],['+','Add']],[['<','lt'],['>','gt'],['<=','lte'],['>=','gte']]];\n      // For asciimath, some fixed translations apply (like pi->Math.PI) etc ..                                          \n        tok=tok.map(t=>(t[0]!=5)?t:[].concat.apply([],syntax).filter(x=>x[0]==t[1]).length?[5,[].concat.apply([],syntax).filter(x=>x[0]==t[1])[0][1]]:t); \n      // Now the token-stream is translated recursively.    \n        function translate(tokens) {\n           // helpers : first token to the left of x that is not of a type in the skip list.\n           var left = (x=ti-1,skip=[0])=>{ while(x>=0&&~skip.indexOf(tokens[x][0])) x--; return x; },\n           // first token to the right of x that is not of a type in the skip list.\n               right= (x=ti+1,skip=[0])=>{ while(x<tokens.length&&~skip.indexOf(tokens[x][0])) x++; return x; },\n           // glue from x to y as new type, optionally replace the substring with sub.    \n               glue = (x,y,tp=5,sub)=>{tokens.splice(x,y-x+1,[tp,...(sub||tokens.slice(x,y+1))])},\n           // match O-C pairs. returns the 'matching bracket' position    \n               match = (O=\"(\",C=\")\")=>{var o=1,x=ti+1; while(o){if(tokens[x][1]==O)o++;if(tokens[x][1]==C)o--; x++;}; return x-1;};\n           // grouping (resolving brackets).    \n           for (var ti=0,t,si;t=tokens[ti];ti++) if (t[1]==\"(\") glue(ti,si=match(),6,[[4,\"(\"],...translate(tokens.slice(ti+1,si)),[4,\")\"]]); \n           // [] dot call and new\n           for (var ti=0,t,si; t=tokens[ti];ti++) {\n             if (t[1]==\"[\") { glue(ti,si=match(\"[\",\"]\"),6,[[4,\"[\"],...translate(tokens.slice(ti+1,si)),[4,\"]\"]]); if (ti)ti--;}    // matching []\n             else if (t[1]==\".\") { glue(left(),right()); ti--; }                                                                   // dot operator\n             else if (t[0]==6 && ti && left()>=0 && tokens[left()][0]>=5 && tokens[left()][1]!=\"return\") { glue(left(),ti--) }     // collate ( and [\n             else if (t[1]=='new') { glue(ti,right()) };                                                                           // collate new keyword\n           }\n           // ++ and --\n           for (var ti=0,t; t=tokens[ti];ti++) if (t[1]==\"++\" || t[1]==\"--\") glue(left(),ti);  \n           // unary - and + are handled seperately from syntax ..\n           for (var ti=0,t,si; t=tokens[ti];ti++) \n             if (t[1]==\"-\" && (left()<0 || (tokens[left()]||[4])[0]==4)) glue(ti,right(),5,[\"Element.Sub(\",tokens[right()],\")\"]);   // unary minus works on all types.\n             else if (t[1]==\"+\" && (tokens[left()]||[0])[0]==4) glue(ti,ti+1);                                                      // unary plus is glued, only on scalars.\n           // now process all operators in the syntax list .. \n           for (var si=0,s; s=syntax[si]; si++) for (var ti=s[0][3]?tokens.length-1:0,t; t=tokens[ti];s[0][3]?ti--:ti++) for (var opi=0,op; op=s[opi]; opi++) if (t[1]==op[0]) {\n             // exception case .. \".Normalized\" and \".Length\" properties are re-routed (so they work on scalars etc ..)\n                    if (op[2]==2) { var arg=tokens[left()]; glue(ti-1,ti,5,[\"Element.\"+op[1],\"(\",arg,\")\"]); }                                    \n             // unary operators (all are to the left)      \n               else if (op[2])    { var arg=tokens[right()]; glue(ti, right(), 5, [\"Element.\"+op[1],\"(\",arg,\")\"]); }\n             // binary operators  \n                             else { var l=left(),r=right(),a1=tokens[l],a2=tokens[r]; glue(l,r,5,[\"Element.\"+op[1],\"(\",a1,\",\",a2,\")\"]); ti--; }\n           }\n           return tokens;\n        }        \n      // Glue all back together and return as bound function.  \n        return eval( ('('+(function f(t){return t.map(t=>t instanceof Array?f(t):typeof t == \"string\"?t:\"\").join('');})(translate(tok))+')') );\n      }\n    }\n    \n    \n  // Matrix-free inverses up to 5D. Should translate this to an inline call for readability.\n  // http://repository.essex.ac.uk/17282/1/TechReport_CES-534.pdf  \n    res.prototype.__defineGetter__('Inverse', function(){  \n      return (tot==0)?new this.constructor.Scalar([1/this[0]]):\n             (tot==1)?this.Involute.Mul(this.constructor.Scalar(1/this.Mul(this.Involute)[0])):\n             (tot==2)?this.Conjugate.Mul(this.constructor.Scalar(1/this.Mul(this.Conjugate)[0])):\n             (tot==3)?this.Reverse.Mul(this.Involute).Mul(this.Conjugate).Mul( this.constructor.Scalar(1/this.Mul(this.Conjugate).Mul(this.Involute).Mul(this.Reverse)[0])):\n             (tot==4)?this.Conjugate.Mul(this.Mul(this.Conjugate).Map(3,4)).Mul( this.constructor.Scalar(1/this.Mul(this.Conjugate).Mul(this.Mul(this.Conjugate).Map(3,4))[0])):\n                      this.Conjugate.Mul(this.Involute).Mul(this.Reverse).Mul(this.Mul(this.Conjugate).Mul(this.Involute).Mul(this.Reverse).Map(1,4)).Mul(this.constructor.Scalar(1/this.Mul(this.Conjugate).Mul(this.Involute).Mul(this.Reverse).Mul(this.Mul(this.Conjugate).Mul(this.Involute).Mul(this.Reverse).Map(1,4))[0]));\n    });\n    \n  // If a function was passed in, translate, call and return its result. Else just return the Algebra.  \n    if (fu instanceof Function) return res.inline(fu)(); else return res;  \n  }\n}));\n\n        function add_graph_to_notebook(Algebra){\n            var output = Algebra({p:4,q:1,r:0,baseType:Float64Array},()=>{\n                // When we get a file, we load and display.\n                var canvas;\n                var h=0, p=0;\n                // convert arrays of floats back to CGA elements.\n                var data = [0, {data:[[0.0, 2.9411635777003005, 16.575365852710323, 17.074564383213925, 286.9669726092519, 287.9669726092519, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.840493862435234e-14, 3.1425630951060477e-13, 3.1425630951060477e-13, -3.053009837860811e-13, -3.053009837860811e-13, 0.0, 5.54876769962813e-14, 5.54876769962813e-14, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, -17.74964232619511, -17.382699467468377, 8.296485688625085, 342.5198591328716, 343.5198591328716, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2.902259007827684e-13, 2.407855032312562e-12, 2.407855032312562e-12, 5.051986944797036e-12, 5.051986944797036e-12, 0.0, -5.147867669797527e-12, -5.147867669797527e-12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 10.059454050071313, -1.6362200299057186, 15.910162268580523, 178.00154759216207, 179.00154759216207, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -6.168350758509988e-14, -9.813946149741565e-13, -9.813946149741565e-13, -1.0730731417359889e-13, -1.0730731417359889e-13, 0.0, -6.211136519948125e-13, -6.211136519948125e-13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 7.764469102554897, -9.9621273984128, 13.879219950839312, 175.58185459525632, 176.58185459525632, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.230267640781833e-13, -1.7075155184811227e-12, -1.7075155184811227e-12, -1.2238709013343548e-12, -1.2238709013343548e-12, 0.0, -9.607878145530505e-13, -9.607878145530505e-13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, -18.926001924439852, 18.66851117686507, 13.088120140897335, 438.5028736135998, 439.5028736135998, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.603519071673545e-13, -2.098705025828352e-12, -2.098705025828352e-12, 2.988905430688904e-12, 2.988905430688904e-12, 0.0, 3.0379717063450904e-12, 3.0379717063450904e-12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 3.749118757761751, 17.593179526282555, 16.2529285917069, 293.3667725554661, 294.3667725554661, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.5334960575735212e-13, 2.4923801919406497e-12, 2.4923801919406497e-12, -2.6989847701310645e-12, -2.6989847701310645e-12, 0.0, 5.716800360241692e-13, 5.716800360241692e-13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 7.5366619809381294, -11.169047416322057, 12.75784736438073, 171.6557816879025, 172.6557816879025, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -4.9054708338539577e-14, -6.258324814873026e-13, -6.258324814873026e-13, -5.510588486599378e-13, -5.510588486599378e-13, 0.0, -3.6959548895816597e-13, -3.6959548895816597e-13, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, -11.882567539962201, -20.58822454912611, 8.14589236015862, 315.21298188514646, 316.21298188514646, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -2.018652604271175e-13, -1.6443726826946864e-12, -1.6443726826946864e-12, -4.152808443590151e-12, -4.152808443590151e-12, 0.0, 2.3930505915756095e-12, 2.3930505915756095e-12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, 10.708656243917957, 10.016465066593248, 16.579794819308823, 244.4472436155328, 245.4472436155328, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 7.913914480366857e-15, 1.312110783020395e-13, 1.312110783020395e-13, -8.226728811190655e-14, -8.226728811190655e-14, 0.0, 8.500625900261265e-14, 8.500625900261265e-14, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [0.0, -23.446517822856446, -15.582746926750287, 6.574214450599988, 417.3907477215741, 418.3907477215741, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.213336050959161e-13, 2.1125160300849616e-12, 2.1125160300849616e-12, 5.004740553510541e-12, 5.004740553510541e-12, 0.0, -7.532953549010588e-12, -7.532953549010588e-12, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]].map(x=>x.length==32?new Element(x):x).map(x=>x.length==32?new Element(x):x)}, 65280, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.5804935388215091, 12.142436157116785, 12.187157184042686, 0.21684567109922143, 0.2232174797363051, 0.11657617864712576, -6.691459470249603, -6.700454387901287, 0.32735717645423734, 0.07008896322114992, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 16711680, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6577975512778201, 12.162639898466134, 12.209889688880983, -0.022776734618162132, -0.01595948405271484, 0.12766228693489323, -7.296059327799348, -7.305030810363966, 0.3582056342222838, 0.07591058293103932, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], 255, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.6234240943222209, 12.347746243120179, 12.393667737724648, 0.2145462422514287, 0.22133394901992573, 0.11863608948987134, -6.861510060675536, -6.869938827696023, 0.3384766851975262, 0.07180594782524109, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]].map(x=>x.length==32?new Element(x):x);\n                data = data.map(x=>x.length==32?new Element(x):x);\n                // add the graph to the page.\n                canvas = this.graph(data,{gl:true,conformal:true,grid:true,scale:0.05,useUnnaturalLineDisplayForPointPairs:true});\n                canvas.options.h = h; canvas.options.p = p;\n                // make it big.\n                canvas.style.width = '100%';\n                canvas.style.height = '50vh';\n                return canvas;\n            });\n            element.append(output);\n\n            var a = document.createElement(\"button\");\n            var t = document.createTextNode(\"💾 Save\");\n            a.appendChild(t);\n            function screenshot(){\n                //output.width = 1920;  output.height = 1080;\n                output.update(output.value);\n                output.toBlob(function(blob) {\n                    var url = URL.createObjectURL(blob);\n                    window.open(url, '_blank');\n                });\n            }\n            a.onclick = screenshot\n            var butnelem = element.append(a);\n        }\n        // requirejs works in sphinx, require works in jupyter\n        (requirejs || require)(['Algebra'],function(Algebra){add_graph_to_notebook(Algebra)});\n        "
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "# Test them several time"
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "noise = 0.1\nnpnts = 10\nntests = 1000\n\ncdiff_array = np.zeros((ntests, 2))\nndiff_array = np.zeros((ntests, 2))\nrdiff_array = np.zeros((ntests, 2))\nr_array = np.zeros(ntests)\n\nfor i in range(ntests):\n    # The true circle\n    true_circle = random_circle()\n    point_list = project_points_to_circle([random_conformal_point() for i in range(npnts)], true_circle)\n    point_list = [up(down(P) + noise * random_euc_mv()) for P in point_list]\n\n    # Circle from plane sphere intersection\n    sphere = fit_sphere(point_list)\n    plane = fit_plane(point_list)\n    plane_circle = meet(sphere,plane).normal()\n\n    # The circle from leos method\n    leo_circle = fit_circle(point_list)\n\n    ct, mt, rt = get_circle_in_euc(true_circle)\n    cp, mp, rp = get_circle_in_euc(plane_circle)\n    cl, ml, rl = get_circle_in_euc(leo_circle)\n\n    center_diff = [((cp-ct)**2)[0], ((cl-ct)**2)[0]]\n    cdiff_array[i, :] = center_diff\n    normal_diff = [((mp-mt)**2)[0], ((ml-mt)**2)[0]]\n    ndiff_array[i, :] = normal_diff\n    radius_diff = [((rp-rt)**2), ((rl-rt)**2)]\n    rdiff_array[i, :] = radius_diff\n    r_array[i] = rt\n\nprint('Comparison rms error')\nprint('for circles of average radius')\nprint(np.mean(r_array))\nprint('l     p')\nprint('centres')\nprint(np.sqrt(np.mean(cdiff_array, axis=0)))\nprint(np.sqrt(np.median(cdiff_array, axis=0)))\nprint('normals')\nprint(np.sqrt(np.mean(ndiff_array, axis=0)))\nprint(np.sqrt(np.median(ndiff_array, axis=0)))\nprint('radii')\nprint(np.sqrt(np.mean(rdiff_array, axis=0)))\nprint(np.sqrt(np.median(rdiff_array, axis=0)))",
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "text": "Comparison rms error\nfor circles of average radius\n18.65239054046665\nl     p\ncentres\n[20.37996876 18.32594386]\n[1.18569829 1.05548757]\nnormals\n[1.44558323 1.38409866]\n[1.94952453 0.39131832]\nradii\n[9.26278359 9.80583646]\n[0.5003353  0.44704945]\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "",
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3",
      "language": "python"
    },
    "language_info": {
      "name": "python",
      "version": "3.7.4",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "pygments_lexer": "ipython3",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    },
    "gist": {
      "id": "",
      "data": {
        "description": "circle_fitting",
        "public": true
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}