What is being done in here ? (Used math recognition)

Question

I know this isn't exactly programming related per se, but programmers are the most probable of all people who will recognize this maybe.

I have the following (X and Y are arrays, both with 3 elements), and I cannot recognize (although it reminds me of a few things, but none quite!) what is being done here. Does it ring any bells for anyone else ?

I gather you can disregard the lower part; the upper should probably give it away ... but I still cannot see it.

At first it reminded me of linear interpolation in 3d space ...

  SUBROUTINE TRII(X,Y,XR,YR)
DIMENSION X(3),Y(3)

D=X(1)*(X(2)**2-X(3)**2)+
 >    X(2)*(X(3)**2-X(1)**2)+
 >    X(3)*(X(1)**2-X(2)**2)

D1=Y(1)*(X(2)*X(3)**2-X(3)*X(2)**2)+
 >     Y(2)*(X(3)*X(1)**2-X(1)*X(3)**2)+
 >     Y(3)*(X(1)*X(2)**2-X(2)*X(1)**2)

D2=Y(1)*(X(2)**2-X(3)**2)+
 >     Y(2)*(X(3)**2-X(1)**2)+
 >     Y(3)*(X(1)**2-X(2)**2)

D3=X(2)*(Y(3)-Y(1))+
 >     X(1)*(Y(2)-Y(3))+
 >     X(3)*(Y(1)-Y(2))

A=D1/D
B=D2/D
C=D3/D

YR=A+B*XR+C*XR**2

RETURN
END

  SUBROUTINE TRIM(X,Y,XR,YR,XM,YM)
DIMENSION X(3),Y(3)

D=X(1)*(X(2)**2-X(3)**2)+
 >    X(2)*(X(3)**2-X(1)**2)+
 >    X(3)*(X(1)**2-X(2)**2)

D1=Y(1)*(X(2)*X(3)**2-X(3)*X(2)**2)+
 >     Y(2)*(X(3)*X(1)**2-X(1)*X(3)**2)+
 >     Y(3)*(X(1)*X(2)**2-X(2)*X(1)**2)

D2=Y(1)*(X(2)**2-X(3)**2)+
 >     Y(2)*(X(3)**2-X(1)**2)+
 >     Y(3)*(X(1)**2-X(2)**2)

D3=X(2)*(Y(3)-Y(1))+
 >     X(1)*(Y(2)-Y(3))+
 >     X(3)*(Y(1)-Y(2))

A=D1/D
B=D2/D
C=D3/D

XR=-B/(2.*C)
YR=A+B*XR+C*XR**2

XM=XR
IF(XR.GT.X(1).OR.XR.LT.X(3))XM=X(1)
YM=A+B*XM+C*XM**2
IF(YM.LT.Y(1))XM=X(1)
IF(YM.LT.Y(1))YM=Y(1)

RETURN
END

">" is a continuation sign.

Answer 1

I'm not sure what language this is, but it's clear that this is some sort of solver for quadratic equations. The XR and YR expressions are a dead giveaway:

XR = -B / (2.*C)
YR = A + B*XR + C*XR**2

Without knowing what the X(1..3) and Y(1..3) expressions are, however, it's not going to be possible to infer too much more about what the A/B/C coefficients represent, however. Lots of things use quadratic equations -- area of a circle given the radius, intensity of light at a given distance, et cetera. More contextual data is required.

---

Update: The OP indicated that he can't be too much more specific for secrecy reasons. Here are some hints, though:

• What does the subroutine return? How are those results used later on? That may lead to better insights.

• It appears that Y(1) is some sort of magic lower bound for the result of this computation. Notice that if YM is less than Y(1), then both XM and YM are set to X(1) and Y(1), respectively.

• The "D" expressions look like this, in more natural syntax:

d = x1 * [x2^2 - x3^2] + x2 * [x3^2 - x1^2] + x3 * [x1^1 - x2^2]
  d1 = y1 * [x2*x3^2 - x3*x2^2] + y2 * [x3*x1^2 - x1*x3^2] + y3 * [x1*x2^2 - x1*x2^2]
  d2 = y1 * [x2^2 - x3^2] + y2 * [x3^2 - x1^2] + y3 * [x1^2 - x2^2]
  d3 = x2 * [y3 - y1] + x1 * [y2 - y3] * x3 * [y1 - y2]

• This looks very much like some sort of matrix operation; D is almost certainly for "determinant". But there are other things that have the same mathematical relationship.

Answer 2

This is a way to solve linear equation systems, specifically cramers rule. Also have a look at the rule of sarrus. After that, you seem to construct a quadratic equation out of it.

Answer 3

D is the determinant of the matrix:

        | x(1) x(1)² 1 |
D = det | x(2) x(2)² 1 |
        | x(3) x(3)² 1 |

In D1, the rightmost column has been replaced with Y:

         | x(1) x(1)² Y(1) |
D1 = det | x(2) x(2)² Y(2) |
         | x(3) x(3)² Y(3) |

In D2, and D3 it's the first and second columns, respectively. Is it easier to recognize now? Looks a lot like using Cramer's rule to solve a linear equation to me.

Edit: To be more precise: (A, B, C) is the solution to the system:

A + x(1)*B + x(1)²*C = Y(1)
A + x(2)*B + x(2)²*C = Y(2)
A + x(3)*B + x(3)²*C = Y(3)

YR is the square of the solution to the quadratic equation (nb, different x!):

C*x² + B*x + A = 0

I feel like this should be obvious now, but I can't quite grasp it...

Answer 4

This code represents a kind of interpolation/quadratic curve fitting on three 2d points together with a way to compute the minimum or maximum value of such a fitted quadratic within the interval itself. I guess that TRII stands for triple (point)-interpolation and TRIM stands for triple (point) minimum or maximum.

To be more precised TRII solves the problem :- find a quadratic curve that passes through the points (x1,y1),(x2,y2) and (x3,y3) in the form Y=A+B*X+C*X^2 and compute the Y value of the quadratic at the point XR and return as YR. This is basically a way to interpolate smoothly between three 2d points. It is often used to find a better approximation for the max or min value of a set of discrete data points.

All the D, D1, D2, D3 stuff is to solve the matrix equation:

(1 X1 X1^2) *(A) = (Y1)

(1 X2 X2^2) *(B) = (Y2)

(1 X3 X3^2) *(C) = (Y3)

using Cramers rule as mentioned in one of the other comments, D is the matrix determinant and D1, D2, D3 are co-factors.

TRIM again computes the quadratic Y=A+B*X+C*X^2 and then finds a max/min of this quadratic (XM, YM). This is done by initially finding the point where the quadratic has a turning point: if F(X)=A+B*X+C*X^2, F'(XR)=B+2*C*XR=0, or XR=-B/2*C, YR=A+B*XR+C*XR^2. There is then some logic to force the returned XM, YM min or max values to lie within certain bounds.

The code:

XM=XR . . . IF(YM.LT.Y(1))YM=Y(1)

Is a little weird since if we assume that GT and LT mean greater than and less than respectively then we need to assume that X3'<'X1 otherwise the condition (XR.GT.X(1).OR.XR.LT.X(3)) is trivial and XM,YM are set to X1, Y1.

So X3'<'X1 and the condition says that if the quadratics max/min value is outside the interval (X1,X3) then set (XM,YM) to (X1, Y1) as before. If not then if Y1 is above the min/max value in Y then again set (XM,YM) to (X1, Y1).

It is hard to understand what this means and I suspect the code may be wrong! Any thoughts?

Ivan

Answer 5

The code run as follows

Routine TRII takes as input the coordinates of three points (x,y) and interpolates a parabola using Lagrange interpolation. Also takes as input the coordinate XR. Returns in YR the value at XR for the interpolating parabola. I guess the name of the routine comes from "TRI" (Croatian for "three" (points)) and "I" for Interpolation.

Routine TRIM also calculates the same parabola, and returns the minimun value of the function in the interval {X(1),X(3)}

(I "really" executed the program) >)

Note that this is FORTRAN code and the parameters are passed by reference, so the results are returned back in the same parameters (very odd!)

EDIT: Just an old anecdote

Many years ago, before the PC era, I worked as a junior programmer in the IBM5110 platform. The machine had two interpreters (Basic and APL) and no compiler. We used APL for calculation intensive tasks, but mostly Basic for commercial applications. The Basic interpreter had very limited control structures (IF - without else, FOR , GOTO/ GOSUB and GOTO/GOSUB .. ON ..).

The GOSUB .. ON had a nasty syntax like this

 GOSUB 1025,2033,3008 ON exp

where 1025,2033 and 3008 were line numbers (the language didn't support labels) and "exp" was any formula or variable evaluated (in this example) 1,2 or 3. Upon execution the program transferred control to line 1025 if exp was evaluated as 1, and so on.

I used it all the time, calculating a Vandermonde polynomial on a variable to get the desired result (1, 2, 3 ...). I was young, and code maintenance was not a big issue for me ... until my boss saw what I was doing. He was very kind, and decided not to kill me.

BTW, the machines had 48KB RAM, no hard disk, and 1MB floppy disks. We developed a lot of commercial apps (flights reservations, general ledger, billing, stocks trading ...)