Moving parabolic
approximation (in C)




This is the functionMovingParabApprox converted into C. The conversion was done by the OMtoC function (written in Maple). Except for wrapping around two long lines, the result shown here is the exact output of the converter.

void MovingParabApprox (int n, int m, float x[], float a0[],
  float a1[], float a2[] )
{
  int i; float sx; float sjx; float sj2x;
  if( m < 3 )
     { printf( "An irrecoverable error has occurred\n" );
       exit(1);
     }
  if( 1 <= n )
     { a0[0] = x[0];
       a1[0] = 0;
       a2[0] = 0;
     }
  if( 2 <= n )
     { a0[1] = x[1];
       a1[1] = x[1]-x[0];
       a2[1] = 0;
     }
  /* The definitions of sx, sjx, and sj2x are:
 m - 1                m - 1                   m - 1
 -----                -----                   -----
  \                    \                       \     2
   )   x[i - j] = sx,   )   j x[i - j] = sjx,   )   j  x[i - j] = sj2x
  /                    /                       /
 -----                -----                   -----
 j = 0                j = 0                   j = 0 */
  if( 3 <= n )
     { sx = x[0]+x[1];
       sjx = x[0];
       sj2x = x[0];
     }
  for( i=3; i <= (m<n?m:n); i += 1 )
    { sj2x += 2*sjx+sx;
      sjx += sx;
      sx += x[i-1];
      { float t8; float t11; float t12; float t13;
          float t14; float t15;
        t15 = i*sjx;
        t8 = i*i;
        t14 = sx*t8;
        t13 = 1./((i+1)*(i+2)*i);
        t12 = 6*sjx-(3*i-2)*sx;
        t11 = t13/((i-2)*(i-1));
        a1[i-1] = ((-180*i+180)*sj2x-(-192*t8+360*i-132)*
          sjx-(-126*t8-36+(126+36*t8)*i)*sx)*t11;
        a0[i-1] = (9*t14-36*t15+30*sj2x+3*t12)*t13;
        a2[i-1] = (30*t14-180*t15+180*sj2x+30*t12)*t11;
      }
    }
  for( i=m+1; i <= n; i += 1 )
    { sj2x -= x[i-m-1]*(m-1)*(m-1);
      sjx -= x[i-m-1]*(m-1);
      sx -= x[i-m-1];
      sj2x += 2*sjx+sx;
      sjx += sx;
      sx += x[i-1];
      { float t8; float t11; float t12; float t13;
          float t14; float t15;
        t15 = m*sjx;
        t8 = m*m;
        t14 = sx*t8;
        t13 = 1./((m+1)*(m+2)*m);
        t12 = 6*sjx-(3*m-2)*sx;
        t11 = t13/((m-2)*(m-1));
        a1[i-1] = ((-180*m+180)*sj2x-(-192*t8+360*m-132)*
          sjx-(-126*t8-36+(126+36*t8)*m)*sx)*t11;
        a0[i-1] = (9*t14-36*t15+30*sj2x+3*t12)*t13;
        a2[i-1] = (30*t14-180*t15+180*sj2x+30*t12)*t11;
      }
    }
  return;
}