Math.hpprgm

// JMB 2018-08
#pragma mode( separator(.,;) integer(h32) )

EXPORT Polynomial_Interpolation(data,x)
// data:[[x1,y1], x:abscissa → interpolated ordinate y
//       [x2,y2]
//       .......
//       [xn,yn]]  
BEGIN
  RETURN POLYEVAL(polynomial_regression(data,rowDim(data)-1),x);
END;

EXPORT Polynomial_Interpolation_2D(data,x,y)
// data:[[0,   x1,  x2, ...  xn], x, y → interpolated z
//       [y1, z11, z12, ... z1n]
//       [y2, z21, z22, ... z2n]
//       .......................
//       [ym, zm1, zm2, ... zmn]]  
BEGIN
  LOCAL c := colDim(data);
  LOCAL r := rowDim(data);
  LOCAL d := [[0,0],[0,0]];
  LOCAL z := [[0,0],[0,0]];
  LOCAL k;
  d(-1) := data({{1,2},{1,c}});
  FOR k FROM 2 TO r DO
    d(-2) := data({{k,2},{k,c}});
    z(k-1,2) := Polynomial_Interpolation(d,x);
  END;
  z(-1) := TRN(data({{2,1},{r,1}}));
  RETURN Polynomial_Interpolation(z,y);
END;

EXPORT Trapezoidal_Rule(data)
// data:[[x1,y1] → Integral from x1 to xn
//       [x2,y2]
//       .......
//       [xn,yn]]  
BEGIN
  LOCAL k;
  LOCAL a := 0;
  FOR k FROM 1 TO rowDim(data)-1 DO
    a := a+(data(k+1,2)+data(k,2))*(data(k+1,1)-data(k,1));
  END;
  RETURN a/2;
END;

EXPORT REC_to_CYL(rec)
// rec:rectangular coordinates [x,y,z] → cylindrical coordinates [ρ,φ,z]
BEGIN
  LOCAL cyl := polar_coordinates(rec(1),rec(2));
  RETURN [cyl(1),cyl(2),rec(3)];
END;

EXPORT CYL_to_REC(cyl)
// cyl:cylindrical coordinates [ρ,φ,z] → rectangular coordinates [x,y,z]
BEGIN
  LOCAL rec := rectangular_coordinates(cyl(1),cyl(2));
  RETURN [rec(1),rec(2),cyl(3)];
END;
  
EXPORT SPH_to_CYL(sph)
// sph:spherical coordinates [r,θ,φ] → cylindrical coordinates [ρ,φ,z]
BEGIN
  LOCAL cyl := rectangular_coordinates(sph(1),sph(2));
  RETURN [cyl(2),sph(3),cyl(1)];
END;
 
EXPORT CYL_to_SPH(cyl)
// cyl:cylindrical coordinates [ρ,φ,z] → spherical coordinates [r,θ,φ]
BEGIN
  LOCAL sph := polar_coordinates(cyl(3),cyl(1));
  RETURN [sph(1),sph(2),cyl(2)];
END; 
  
EXPORT REC_to_SPH(rec)
// rec:rectangular coordinates [x,y,z] → spherical coordinates [r,θ,φ]
BEGIN
  RETURN CYL_to_SPH(REC_to_CYL(rec));
END;

EXPORT SPH_to_REC(sph)
// sph:spherical coordinates [r,θ,φ] → rectangular coordinates [x,y,z]
BEGIN
  RETURN CYL_to_REC(SPH_to_CYL(sph));
END;

EXPORT Optimal_Scale(Minimum, Maximum, NumDivTry)
// Minimum, Maximum, NumDivTry (list) → {Bottom, Top, Step, Number of divisions}
BEGIN
  LOCAL BottomTry, TopTry, StpTry, EfficiencyTry; 
  LOCAL Bottom, Top, Stp, Efficiency, NumDiv, I, K;
  Efficiency := 0;
  FOR K FROM 1 TO SIZE(NumDivTry) DO
    StpTry := 10^FLOOR(LOG((Maximum-Minimum)/NumDivTry(K)));
    I := 1;
    REPEAT
      BottomTry := StpTry*FLOOR(Minimum/StpTry);
      TopTry := BottomTry+StpTry*NumDivTry(K);
      IF Maximum ≤ TopTry THEN
        BREAK;
      END;
      I := I+1;
      IF (I MOD 4) == 0 THEN
        StpTry := StpTry*1.25;
      ELSE
        StpTry := StpTry*2;
      END;
    UNTIL false;
    EfficiencyTry := (Maximum-Minimum)/(TopTry-BottomTry);
    IF EfficiencyTry > Efficiency THEN
      Efficiency := EfficiencyTry;
      Bottom := BottomTry;
      Top := TopTry;
      Stp := StpTry;
      NumDiv := NumDivTry(K);
    END; 
  END;  
  RETURN {Bottom,Top,Stp,NumDiv};
END;

EXPORT Decimal_to_Roman(d)
// d:decimal number (from 1 to 3999) → roman number (string)
BEGIN
  LOCAL rl := {"M","CM","D","CD","C","XC","L","XL","X","IX","V","IV","I"};
  LOCAL dl := {1000,900,500,400,100,90,50,40,10,9,5,4,1};
  LOCAL k,rn,dn;
  IF 0 < d < 4000 THEN
    rn := "";
    dn := d;
    FOR k FROM 1 TO 13 DO
      WHILE dn ≥ dl(k) DO 
        rn := rn+rl(k);
        dn := dn-dl(k);
      END;
    END;
    RETURN rn;
  ELSE
    RETURN "Error:Invalid number";
  END;
END;

EXPORT Roman_to_Decimal(r)
// r:roman number (string) → decimal number
BEGIN
  LOCAL dl := {1,5,10,50,100,500,1000};
  LOCAL rl, k, rn, dn, cur, prev, sz;
  IF TYPE(r) ≠ 2 THEN
    RETURN "Error:String expected";
  ELSE
    rn := UPPER(r);
    sz := SIZE(rn);
    rl := MAKELIST(POS(ASC("IVXLCDM"),rn(k)),k,1,sz);
    IF NOT ΠLIST(rl) THEN
      RETURN "Error:Invalid number"; // non-roman characters
    ELSE
      dn := 0;
      FOR k FROM 1 TO sz DO
        cur := dl(rl(k));
        dn := dn+IFTE(prev < cur, -prev, prev);
        prev := cur;
      END;
      RETURN dn+cur;
    END;
  END;
END;