Mercury Has Landed: Arbitrary-Precision Math for JWCEssentials

I am happy to announce that Mercury has now joined the JWCEssentials family.

Mercury is a native C library for arbitrary-precision floating-point arithmetic. It is designed around a simple but powerful idea: instead of forcing large numbers through decimal-first thinking, Mercury stores them in a machine-native base-2³² positional format.

That means each “digit” of the number is a 32-bit unsigned integer. On modern hardware, this is a very natural way to build large-number arithmetic, because two 32-bit words can be multiplied with room for the result and carry inside a 64-bit operation. It is not base 10 dressed up for the machine. It is the machine’s own kind of arithmetic brought closer to the surface.

Why Mercury?

One of the most interesting lessons from writing big-number software is that representation matters. Decimal is excellent for people. Binary and hexadecimal are excellent for machines. Mercury leans into that distinction.

A Mercury number is stored as a sign flag, an exponent, and a sequence of 32-bit mantissa words. The number can then be viewed cleanly through hexadecimal, where each 32-bit word maps perfectly to eight hex digits. This gives hexadecimal a special role: it becomes a transparent window into the native representation, rather than a lossy-looking conversion layer.

That makes Mercury especially useful for the kinds of topics I want to write more about: arbitrary precision, numerical representation, computation, multiplication, division, square roots, powers, logs, and eventually the bridge between numerical engines and computer algebra systems.

The Native Core and the Managed Wrapper

Mercury itself lives as a native C library under:

JWCEssentials/C/Mercury

The managed wrapper lives in:

JWCEssentials/Project/JWCEssentials.net

The wrapper assembly is named Mercury.net. I chose that name intentionally. Mercury should be able to stand on its own as a numerical engine without being over-branded. JWCEssentials can host it, document it, and build around it, but the library itself deserves a clean identity.

On the .NET side, the wrapper exposes the native engine through UltraNumber, giving C# code a much friendlier surface: operator overloads, parsing, string conversion, precision management, and thread-local stack handling.

The Scratch Stack Model

One of Mercury’s design choices that I especially like is that arithmetic functions do not casually allocate temporary memory from the heap. Instead, Mercury uses a caller-supplied scratch stack.

That design keeps temporary allocation explicit, deterministic, and fast. A function asks the stack for temporary space, performs the calculation, and then releases that space in reverse order. This makes the native engine easier to reason about and keeps the door open for GPU-minded execution models, where memory discipline matters.

In other words, Mercury is not just a big-number library. It is also an experiment in making the “nuts and bolts” of arbitrary-precision computation visible and teachable.

A Base-2³² View of Math

Most of us are trained to think in base 10. That is natural, cultural, and useful. But the machine does not owe base 10 any special loyalty.

When Mercury prints a number in hexadecimal, it is not merely choosing a programmer-friendly display format. It is showing the number in a form that lines up directly with the internal representation. Eight hex digits correspond exactly to one 32-bit word.

That makes hexadecimal a kind of inspection window. It lets us look at a large computed value without immediately paying the cost, complexity, and conceptual distraction of decimal conversion.

For example, π begins in hexadecimal as:

3.243F6A8885A308D313198A2E03707344A4093822...

That may look unusual if decimal is the only familiar lens, but it is not less mathematical. It is simply a different coordinate system for the same value.

Where This Is Going

Mercury gives me a stable base for a series of articles I have wanted to write for a long time.

I want to talk about multiplication at the limb level. I want to show how large-number division works. I want to explore square roots, powers, logarithms, and constants such as π and e. I also want to connect those topics to the larger question of computer algebra systems: where symbolic mathematics ends, where numerical computation begins, and how the two can support one another.

I am also looking forward to making some of this visual. CrystalCatalyst gives me a way to draw the behavior of algorithms, not just describe them. A geometric sine and cosine solver, arbitrary-precision constants, and native hexadecimal output all belong to the same larger story: math is not only something we calculate. It is something we can watch unfold.

Status

Mercury is now building on both Ubuntu and Windows, with the native library and managed wrapper in place. The project is young, but the foundation is real.

This is the beginning of a new numerical-computation thread inside JWCEssentials. Mercury is the native body. Mercury.net is the managed bridge. UltraNumber is the friendly C# face.

And for the blog, this opens the door to something I am very excited about: showing big-number math not as a black box, but as a living set of understandable operations.

Mercury has landed.

JWCEssentials on GitHub

The Arbitrary Nature of Base 10 and the Elegance of Base $2^{32}$

We are universally conditioned to think in base 10, but mathematically speaking, all number systems are created equal. A base is simply a human (or mechanical) convention for representing a value. Integers represent exact values regardless of the base they are constructed in—there is never “room left over” or “unused space” in the absolute value itself.

However, when we map math onto silicon, the base we choose matters immensely.

Under the hood, most modern computers chew through data using 32-bit integers. This architecture is virtually flawless for building custom, arbitrary-precision math libraries. Why? Because all fundamental mathematical primitives can be constructed safely using 64-bit hardware operations. You can multiply two 32-bit “digits” together, and the resulting value—along with its carry—fits perfectly into a 64-bit register. No overflow. No lost data. For a machine, base $2^{32}$ is the ideal foundation for a number system.

The Problem with Decimal Conversion

When dealing with massive, precise fractional numbers natively in binary, converting the output to base 10 is an expensive and ugly process. You lose the exact 1-to-1 mapping of the machine’s memory, trading computational elegance for human readability.

The Hexadecimal Window

Hexadecimal (base 16) solves this. Because 16 divides flawlessly into 32 (eight hex digits per 32-bit integer), Hex provides a perfect, transparent window into the native bit-wise estimation of the machine. It offers 100% coverage of the number with zero alignment loss.

20260627_162840: Stay tuned for a math library with the entire mathematical set (plus,minus,times,divide,power,log,root) implemented algorithmically and the continuance of Truth in the Flip

To demonstrate what native base $2^{32}$ mathematics looks like when viewed through a Hexadecimal lens, I let the engine stretch its legs. Here is an absurdly precise calculation of Pi, untouched by base-10 conversion algorithms (truncated for space):

3.243F6A8885A308D313198A2E03707344A4093822299F31D0082EFA98EC4E6C89452821E638D01377BE5466CF34E90C6CC0AC29B7C97C50DD3F84D5B5B54709179216D5D98979FB1BD1310BA698DFB5AC2FFD72DBD01ADFB7B8E1AFED6A267E96BA7C9045F12C7F9924A19947B3916CF70801F2E2858EFC16636920D871574E69A458FEA3F4933D7E0D95748F728EB658718BCD5882154AEE7B54A41DC25A59B59C30D5392AF26013C5D1B023286085F0CA417918B8DB38EF8E79DCB0603A180E6C9E0E8BB01E8A3ED71577C1BD314B2778AF2FDA55605C60E65525F3AA55AB945748986263E8144055CA396A2AAB10B6B4CC5C341141E8CEA15486AF7C72E993B3EE1411636FBC2A2BA9C55D741831F6CE5C3E169B87931EAFD6BA336C24CF5C7A325381289586773B8F48986B4BB9AFC4BFE81B6628219361D809CCFB21A991487CAC605DEC8032EF845D5DE98575B1DC262302EB651B8823893E81D396ACC50F6D6FF383F442392E0B4482A484200469C8F04A9E1F9B5E21C66842F6E96C9A670C9C61ABD388F06A51A0D2D8542F68960FA728AB5133A36EEF0B6C137A3BE4BA3BF0507EFB2A98A1F1651D39AF017666CA593E82430E888CEE8619456F9FB47D84A5C33B8B5EBEE06F75D885C12073401A449F56C16AA64ED3AA62363F77061BFEDF72429B023D37D0D724D00A1248DB0FEAD349F1C09B075372C980991B7B25D479D8F6E8DEF7E3FE501AB6794C3B976CE0BD04C006BAC1A94FB6409F60C45E5C9EC2196A246368FB6FAF3E6C53B51339B2EB3B52EC6F6DFC511F9B30952CCC814544AF5EBD09BEE3D004DE334AFD660F2807192E4BB3C0CBA85745C8740FD20B5F39B9D3FBDB5579C0BD1A60320AD6A100C6402C7279679F25FEFB1FA3CC8EA5E9F8DB3222F83C7516DFFD616B152F501EC8AD0552AB323DB5FAFD23876053317B483E00DF829E5C57BBCA6F8CA01A87562EDF1769DBD542A8F6287EFFC3AC6732C68C4F5573695B27B0BBCA58C8E1FFA35DB8F011A010FA3D98FD2183B84AFCB56C2DD1D35B9A53E479B6F84565D28E49BC4BFB9790E1DDF2DAA4CB7E3362FB1341CEE4C6E8EF20CA

N’th-Dimensional Interpolation Revisited: When a Point Becomes a Sample

When I first wrote about N’th-dimensional interpolation on an array, the goal was straightforward:
given a coordinate with fractional parts, find the neighboring integer coordinates and interpolate
between them. In two dimensions this becomes bilinear interpolation. In three dimensions it becomes
trilinear interpolation. In N dimensions, the same idea generalizes naturally.

The core pattern is simple:

coordinate
    → integer base coordinate
    → fractional offset per dimension
    → collect 2^n neighboring corner values
    → fold those values through interpolation
    → final interpolated value

That original version worked by treating each corner as a point. For each generated corner coordinate,
the array was sampled directly:

value = inputArray[cornerCoordinate]

In C# terms, the key line was essentially:

flat[i] = inputArray.GetValue(interpCoords);

That line was correct, but it also hid something important. It made a quiet assumption:
a coordinate sample means one array cell.

The new realization is that this does not have to be true.

The Sampling Seam

The important change is replacing direct array access with a sample delegate:

flat[i] = sample(inputArray, interpCoords);

This small change opens the algorithm up. Each corner no longer has to mean “read one value from
the array.” Each corner can now mean “sample this location according to some rule.”

The original behavior still exists as the default sample:

SampleDefault(array, coords):
    return array[coords]

But once sampling is abstracted, the interpolation algorithm becomes more than a point interpolator.
It becomes a framework where the meaning of a sample can be changed.

PointSample     → read one cell
BoxSample       → average a local rectangle
CubeSample      → average a local cube
HyperBoxSample  → average a local N-dimensional region
KernelSample    → use weighted neighborhood logic

This is the conceptual upgrade:

Interpolation does not have to interpolate points.
It can interpolate samples.

The Half-Pixel Thought

The discovery that led me back to this code was the idea of a half-pixel.
In two-dimensional image terms, a half-pixel location can be thought of as the space between neighboring
pixels. At exactly halfway between four pixels, ordinary bilinear interpolation naturally averages the
surrounding 2×2 rectangle.

That gives a helpful way to think about the coordinate:

(x, y)       → sample at the pixel/cell position
(x+.5, y+.5) → sample halfway into the neighboring rectangle

In N dimensions, the same idea generalizes:

coords[d] + 0.5

This shifts the sampling coordinate by half a cell in each dimension. Then the existing interpolation
logic does what it already knows how to do: it finds the surrounding 2^n corners and folds them down
into a final value.

In 2D, this means the half-shift samples across a rectangle.
In 3D, it samples across a cube.
In N dimensions, it samples across a hyper-rectangle.

Point Sample vs. Region Sample

There are now two related but distinct ideas:

1. Shift the coordinate
   coords → coords + 0.5

2. Change the sample meaning
   point sample → region/kernel sample

The coordinate shift changes where interpolation happens.
The sample delegate changes what each corner means.

Together, they create a very flexible structure:

coordinate transform
    → corner generation
    → sample delegate
    → interpolation fold

That separation matters. The interpolator does not need to know whether a sample is a point,
a rectangle, a cube, or a weighted neighborhood. It only needs a value for each corner.
The sample function owns the meaning of that value.

Pseudo-Code: The Sample Delegate

The delegate idea can be expressed like this:

SampleDelegate(array, coords):
    return some value of type T from the array at or around coords

The original point sample:

PointSample(array, coords):
    return array[coords]

A simple 2D box sample might look like:

BoxSample2D(array, coords):
    sum = 0
    count = 0

    for dy in 0..1:
        for dx in 0..1:
            p = clamp(coords + (dx, dy))
            sum += array[p]
            count += 1

    return sum / count

A 3D cube sample follows the same pattern:

CubeSample3D(array, coords):
    sum = 0
    count = 0

    for dz in 0..1:
        for dy in 0..1:
            for dx in 0..1:
                p = clamp(coords + (dx, dy, dz))
                sum += array[p]
                count += 1

    return sum / count

And the N-dimensional version is the natural continuation:

HyperBoxSampleND(array, coords):
    sum = 0
    count = 0

    for each offset in all binary offsets for N dimensions:
        p = clamp(coords + offset)
        sum += array[p]
        count += 1

    return sum / count

For N dimensions, the number of offsets in a 2-wide hyper-box is:

2^n

That mirrors the interpolation itself, which also gathers 2^n neighboring corner values.
This symmetry is part of what makes the idea feel natural.

The Updated Interpolation Shape

With the sampling delegate in place, the high-level interpolation algorithm becomes:

Interpolate(array, coords, interpolator, sample):
    split coords into base coordinates and fractional q values

    for each corner among 2^n corners:
        cornerCoord = baseCoord + cornerOffset
        flat[i] = sample(array, cornerCoord)

    while more than one value remains:
        fold values together using q for the current dimension

    return final folded value

The old version is still available by passing the default point sample.
The new version allows richer sampling without rewriting the interpolation fold.

Place for Updated Code

Below is the updated C# implementation.

    public class Nth
    {
        public delegate T SampleDelegate<T>(System.Array inputArray, int[] coords);
        public static T SampleDefault<T>(Array inputArray, int[] coords)
        {
            return (T)inputArray.GetValue(coords);
        }
        
        public delegate T InterpolateDelegate<T>(T a, T b, double q);
        public static double InterpolateDouble(double a, double b, double q)
        {
            return a + (b - a) * q;
        }

        //if you like param arrays here is a nice convenience wrapper
        public static T Interpolate<T>(System.Array inputArray, int[] coords, InterpolateDelegate<T> interpol,
            bool half = false, SampleDelegate<T>? sample = null)
        {
            double[] newCoords = new double[coords.Length];
            for (int i=0; i<coords.Length; i++)
            {
                newCoords[i] = coords[i];
            }
            
            return Interpolate(inputArray, newCoords, interpol, half, sample);
        }

        public static T Interpolate_HalfUnit<T>(System.Array inputArray, int[] coords, InterpolateDelegate<T> interpol,
            SampleDelegate<T>? sample = null)
        {
            return Interpolate(inputArray, coords, interpol, true, sample);
        }
        
        public static T Interpolate<T>(System.Array inputArray, double[] coords, InterpolateDelegate<T> interpol, bool half = false, SampleDelegate<T>? sample = null)
        {
            int dimension;
            int numDimensions = coords.Length;

            if (inputArray.Rank != numDimensions)
                throw new System.ArgumentException("inputArray and coords must have the same number of dimensions");
            
            if (sample == null) 
                sample = SampleDefault<T>;
            
            int stackHeight = 1 << numDimensions;
            
            T[] flat = new T[stackHeight];

            int[] baseCoords = new int[numDimensions];
            int[] interpCoords = new int[numDimensions];
            
            double[] _q = new double[numDimensions];
            if (!half)
            {
                for (dimension = 0; dimension < numDimensions; dimension++)
                {
                    baseCoords[dimension] = (int)Math.Floor(coords[dimension]);
                    _q[dimension] = coords[dimension] - baseCoords[dimension];
                }
            }
            else
            {
                for (dimension = 0; dimension < numDimensions; dimension++)
                {
                    double shifted = coords[dimension] + 0.5;
                    if (shifted >= inputArray.GetLength(dimension))
                        shifted = inputArray.GetLength(dimension) - 1;

                    baseCoords[dimension] = (int) Math.Floor(shifted);
                    _q[dimension] = shifted - baseCoords[dimension];
                }
            }
            
            for (int i = 0; i < stackHeight; i++)
            {
                int ii = i;
                
                for (dimension = 0; dimension < numDimensions; dimension++)
                {
                    int p =  baseCoords[dimension] + (ii % 2);
                    if (p >= inputArray.GetLength(dimension)) p = inputArray.GetLength(dimension) - 1;

                    interpCoords[dimension] = p;

                    ii >>= 1;
                }

                flat[i] = (T) sample(inputArray, interpCoords);
            }

            int foldedStackHeight = stackHeight;
            int dim = numDimensions-1;

            while (foldedStackHeight != 1)
            {
                foldedStackHeight >>= 1;
                for (int position = 0; position < foldedStackHeight; position++)
                {
                    flat[position] = interpol(flat[position], flat[position + foldedStackHeight], _q[dim]);
                    flat[position + foldedStackHeight] = default(T);
                }

                dim--;
            }

            return flat[0];
        }
    }

Why This Matters

The original algorithm answered the question:

How do I interpolate between neighboring points in an N-dimensional array?

The revised version asks a broader question:

What should a sample mean before interpolation happens?

That is a much more powerful question.

For image data, a sample might mean a pixel, a half-pixel blend, or a small rectangle.
For volume data, it might mean a voxel or a cube of voxels.
For procedural fields, it might mean a local kernel.
For generalized numerical arrays, it might mean a neighborhood summary.

The algorithm did not need to become complicated to support this.
It only needed one seam:

array.GetValue(coords)
    → sample(array, coords)

That is the moment where a point becomes a sample.

Closing Thought

This update is exciting to me because it shows the original N’th-dimensional interpolation routine
becoming more general without losing its original simplicity.

The interpolation fold still does the same elegant work:
it reduces 2^n corner values down to one final value.

But now those corner values can carry more meaning.
They can be raw points, half-shifted blends, local regions, or eventually weighted kernels.

In short:

Point → Sample → Region → Kernel

That is a small change in code, but a large change in what the algorithm can express.

JWCEssentials, JWCCommandSpawn, and CrystalCatalystLibrary

I have been hard at work getting several of my foundational libraries ready to share more publicly.
The work has been spread across a few different layers, but the pieces are starting to fit together
in a way that feels important.

At the utility layer, JWCEssentials continues to grow into a practical collection of
building blocks for the larger NewAge ecosystem. One of the pieces I am especially happy with is
JWCCommandSpawn, a command-spawning utility designed to make process execution,
standard input, standard output, and standard error easier to coordinate from managed code.

Command execution sounds simple until you need it to be reliable, reusable, and flexible. Once pipes,
interactive input, output capture, environment setup, and cross-platform behavior enter the picture,
it becomes clear that this is a real infrastructure problem. JWCCommandSpawn is my answer to that
problem: a cleaner way to launch and communicate with external processes without rewriting the same
process-handling logic everywhere.

CrystalCatalystLibrary

The larger milestone is CrystalCatalystLibrary. This project is becoming my
cross-platform native and managed windowing/rendering substrate. It brings together native C++,
managed .NET, pixel presentation, window control, icons, cursors, drag/drop, clipboard work, and
Skia-based rendering.

Recently, CrystalCatalystLibrary gained support for custom cursors, standard cursors, window icons,
window titles, sizing, positioning, and SVG-to-pixel rendering through Skia. The native side now has
platform-specific support on Windows and Linux/X11, while the managed side can work with rendered
pixel data through a common PixData bridge.

One of the most exciting parts is that SVG assets can now become real window assets. An SVG can be
rendered through Skia, converted into PixData, and then used as a window icon, a custom cursor, or
displayed directly in a CrystalCatalyst window.

SVG
  -> Skia renderer
  -> PixData
  -> Window image, icon, or cursor

That means the same visual language can move through the whole stack. A project icon, a cursor, and
a demo animation can all share the same rendering path. That is exactly the kind of seam I like:
one layer of intent flowing through several layers of implementation.

The Managed and Native Bridge

A major part of the design is the bridge between managed and native code. The native layer provides
the platform-specific window behavior. The managed layer gives .NET projects a cleaner surface to use.
Between them is a C-style exported API and a PixData structure for moving pixel buffers across the
boundary.

This gives CrystalCatalystLibrary a practical shape:

Native platform window
  + exported C API
  + managed .NET wrapper
  + PixData interop
  + Skia/SVG rendering
  = CrystalCatalyst surface

This is still evolving, but it is now far enough along that it feels like a real foundation rather
than just an experiment.

Why This Matters

For a long time, I relied on Cairo-style drawing contexts in older projects. Moving toward Skia gives
me a modern drawing foundation that can be shared across applications and libraries. That does create
some version lock, but it also creates consistency. If the whole ecosystem draws through the same
rendering substrate, then icons, controls, cursors, animations, and future UI pieces can all speak
the same graphical language.

JWCEssentials helps with the practical utility layer. JWCCommandSpawn helps with process control.
CrystalCatalystLibrary helps with native windowing, rendering, and presentation. Together, they are
part of the same larger direction: building reusable infrastructure for NewAge, Sigmas, and the
broader set of tools I want to keep developing.

Repositories

The repositories are available from my GitHub profile:

https://github.com/johnwaynecornell?tab=repositories

This is one of those moments where the work starts to feel less like separate pieces and more like a
foundation. The command layer, the native layer, the managed layer, and the rendering layer are all
beginning to line up.

Utilities
  -> Process control
  -> Native windowing
  -> Managed rendering
  -> Reusable application substrate

That is the kind of progress I like to see.