Linear Interpolation Made Convenient in Maxima

load(interpol); pts: [[0,5],[1,1.5],[4.4,3.14159],[10,10]]; f(x) := ''(linearinterpol(pts,varname='x)); g(x) := ''(cspline(pts,varname='x)); wxplot2d([f(t)],[t,-5,15]); wxplot2d([g(t)],[t,-5,15]); wxplot2d([f(t)],[t,0,10]); wxplot2d([g(t)],[t,0,10]);

(I usually put the load statement in a separate cell in wxMaxima so that I am not re-executing the load every time I correct my main code and rerun the cell.  Minor point.)

Maxima can turn a set of points into a function for linear interpolation.  For that matter, it can make a cubic spline out of the same points.  I'm more interested in the linear interpolated function today—the other is for your reference.  There are three main points of interest in the above code:

1. linearinterpol() takes points, sorts them by the first ordinate, and produces an expression.
2. linearinterpol() produces expressions which use the characteristic function, rather than taking cases to make it transparent to other Maxima routines.
3. The use of the double single quote operator ('') and the single quote operator (') helps me to define the function f(x) such that x is a formal variable and not embedded as an expression.
4. linearinterpol() produces an expression which extrapolates beyond the specified points (from negative infinity to positive infinity) and it is left to the user's discretion to determine the limits to which this extrapolation is reasonable or applicable.

In a Nutshell

All you need is a list of 2D points.  There is more than one way to specify these, but I find it easiest to use a list of lists.  The main list elements are all lists of 2 representing x and y ordinates.  The function sorts these points according to the x-coordinate before proceeding.  You can specify the variable which will be used in the expression.  Provided that you precede that variable with a single quote operator, you can use a variable that is elsewhere assigned a value.  More on that below.

Transparent Output

The characteristic function, which is often assigned the Greek letter chi ("χ"), returns either 0 or 1.  It is a function of two variables, a set and a value.  If the value is an element of the set, it returns 1, otherwise it returns 0.  In the context of our linear interpolation expression, we have a linear expression for each interval in our set of points and each of these is multiplied by a characteristic function with a different set.  Since the sets used for the characteristic functions form a partition of the real numbers, we won't have more than one of the terms of our overall expression from linearinterpol() evaluate to a non-zero value at a time.  This expression can be integrated, differentiated, or combined with other expressions without having to take special cases.  Take the monkey off your back, let Maxima do it.

Cool, Useful, but Initially Abstruse Operators

The single quote operator (beside the key) is a way to tell Maxima, "Don't evaluate this thing."  If I write x and x is assigned value, then Maxima will do what is normally very helpful to me, namely, make the substitution for me.  If I don't want Maxima to do that, because I actually want the expression x, not the value that x stands for (because it has been assigned a value), I need to write 'x instead. In contrast, the double single quote operator ('') tells Maxima to replace the given expression with the value of the expression. This way, when I use the delayed assignment operator (:=) it does not delay the execution of the linearinterpol() function, only the evaluation of the expression as a function of x is delayed.  You may recall from my previous post that parameters in function declarations are formal.  Even if x is elsewhere assigned a value, f(x) := some_expression_including_x + x, will still work intuitively.

Here is some sample input followed by the output.  This is the best way to clarify in your mind the use of these operators:

x:2;
y:2*x;
z: 2*'x;
d(x) := z;
d(2);
d(3);
h(x) := ''(z);
h(2);
h(3);

(%o46) 2
(%o47) 4
(%o48) 2*x
(%o49) d(x):=z
(%o50) 2*x
(%o51) 2*x
(%o52) h(x):=2*x
(%o53) 4
(%o54) 6

Extrapolation

The expression that linearinterpol() produces does not protect you from yourself.  If will produce a result from negative to positive infinity without telling you that you have gone beyond the boundaries of your initial data set.  Actually, this is as it should be, in my opinion.  The program gives you the freedom to ask all the what-ifs your little heart desires, but it is up to you to discern what is really true to life and applicable to the particular problem you are trying to solve—whether the extrapolation is truly justified by the data is rightly left as a matter for user judgement.

Making a Double Out of Parts in VB.NET

Recently I picked up on a piece of code I had started several months ago and didn't have time to finish.  I was a bit stumped actually, though I've got it pretty well beat into shape now (maybe more on that later).  I was making a VB.NET type Fraction using BigIntegers as the base type.  As such, the class would theoretically be able to represent any rational number.  But what if I wanted to convert that result to a double?  So, I started into the code necessary to do the job.  I began with the easy part, the cases where the numerator and denominator can both be converted to doubles.  In that case, .NET already knows how to do the needed conversions.  But if the numerator or denominator is too large to be converted to a double, the result of division (mathematically) may still be within the range of a double.  For example,

4.445 × 10⁴⁰⁰ / 4.222 × 10³⁹⁹ = 10.528185693983894

An important finishing piece to solving this programming problem is the following function:

Public Shared Function DoubleByParts(positive As Boolean, exponent As Integer, mantissa As Int64) As Double Dim res As Int64 = 0 If Not positive Then res = &H8000000000000000 Dim biasedExp As Int64 = (exponent + 1023) If biasedExp > &H7FF Then Throw New OverflowException("Exponent too large for conversion to type Double.") If biasedExp < 0 Then Throw New OverflowException("Exponent too small for conversion to type Double.") If (mantissa And &HFFF0000000000000) <> 0 Then Throw New OverflowException("Significand too large for conversion to type Double.") biasedExp <<= 52 res = res Or biasedExp res = res Or mantissa Return BitConverter.Int64BitsToDouble(res) End Function

First, for a good summary of the double format, see Double-precision floating-point format (Wikipedia).

Code Walkthrough

The double type has 64 bits as does the Int64 type.  We begin with initializing our result variable to 0.  Depending on whether the number is indicated as positive or not, we set the sign bit.  The sign bit (bit 63, the most significant bit) needs to be set (1) to indicate a negative number.  For simplicity of use, I wanted to let the caller work in terms of the exponent (base 2) rather than with a biased result.  The 1023 "excess" or bias is one of those values you don't want to have to remember.  This function obviates the need to be overly conscious of that value.

After making our biased exponent (biasedExp), we check that it is in range.  We have exactly 11 bits to work with.  If our result is greater than 11 bits (0x7FF in hex, 11111111111 in binary), then we complain to the caller with an OverflowException.  If the result is negative, then it still does not fit into the 11 bits and we complain with an OverflowException in this case as well.  The mantissa is required to be no more than 52 bits.  That's 13 "nibbles" (4 bits, a single hex digit) or 6 1/2 bytes.  If the mantissa has any bits set in the upper 12 bits, we complain.  By calling the variable mantissa I have already declared my draconian intent with this function.  I will not coddle the caller by shifting things that are in the wrong place to the right place.  If the caller sends a mantissa which is too large, they may have done the shifting wrong or they may have incorrectly included the implicit 1 at the beginning of the "binary scientific notation."  Either way, the caller has made an error and we cannot know which one it is.  If we could know which error they made, it might make sense to correct it silently and move on.  However, a wrong assumption here could be somebody's program's undoing.  Therefore, we just complain and let the caller figure out what to do about it.

The biased exponent is now shifted into place—this, like the 1023 excess, falls into the category of "painful detail we want to encapsulate so we can forget about it."  (This is not ambiguous like the too-large-mantissa.)  In our final steps we "Or" the biased exponent and mantissa into place and use a call to BitConverter.Int64BitsToDouble() to ensure we have the proper format.

Room For Improvement

To make this function more robust, an exception class could be defined, either several different ones (one for each error) or a class with properties and enumerations to allow the caller to efficiently test for the particular condition that caused the "overflow".  I note also that I have used OverflowException where an underflow has occurred.  In program debugging, an accurate error message ("over" versus "under") can provide a clue to the person debugging the program.