# To unbundle, sh this file
echo complextest.cc 1>&2
cat >complextest.cc <<'End of complextest.cc'
#include <stdlib.h>
#include <stdio.h>

#include "mxyzptlk.rsc"

main( int argc, char** argv ) 
{
static double r[] = { 0.0, 0.0 };
int i, count;
DASetup( 2, 6, r );
DA x, y;

if ( argc == 1 ) count = 4;
else             count = atoi( argv[1] );

x.setVariable(0);
y.setVariable(1);

ComplexDA z ( x, y );
ComplexDA w;

w = z;
for ( i = 1; i < count; i++ ){
  w = w*z;
}
w.peekAt();
}
End of complextest.cc
echo concattest.cc 1>&2
cat >concattest.cc <<'End of concattest.cc'
#include <stdlib.h>
#include <stdio.h>
#include "mxyzptlk.rsc"

main() {
DA x, y, u, v, w, z;
LieOperator a, b, c;
static int index[] = { 3, 2 };
double   answer[2];

DASetup( 2, 7, 2 );              

x.setVariable(0);  y.setVariable(1);

u = x*y*y + exp( x + y );  v = cos( y*x*x ) / ( x + 2.0 );  
w = sin(u) * cos(v);       z = exp( u*u*u ) / ( u*v );

a.setVariable( u, 0 ); a.setVariable( v, 1 );

DAFixReferenceAtEnd( a );        
b.fixReference();                
x.setVariable(0);  y.setVariable(1);

u = sin(x) * cos(y);       v = exp( x*x*x ) / ( x*y );

b.setVariable( u, 0 ); b.setVariable( v, 1 );

c = b*a;                         

c.derivative( index, answer );
printf( "Zero: %lf  %lf\n", answer[0], w.derivative( index ) );
printf( "One : %lf  %lf\n", answer[1], z.derivative( index ) );
}
End of concattest.cc
echo coordTest.cc 1>&2
cat >coordTest.cc <<'End of coordTest.cc'
#include "mxyzptlk.rsc"

//***************************************************************
//***************************************************************
//***************************************************************
// 
//
//	Test program.
//

main() {

DASetup( 3, 7 );

coord x( 0.5 );
coord y( 0.3 );
coord z( 0.1 );
DA h, u, v;

static double r[]  = { 1.0,  0.0,  0.0 };
static double s[]  = { 0.33, 0.33, 0.33 };
static int index[] = { 2, 1, 3 };

// First test ....
h = 1.0 / ( cos( exp( x*x*y - 3.0*x*z*z ) ) );
printf( "First test\nThe derivative is: %lf, = %lf\n\n",  
        h.derivative( index ),
        (1.0 / ( cos( exp( x*x*y - 3.0*x*z*z ) ) )).derivative( index ) 
      );

// Second test ....
u = exp( x );
v = exp( x + y + z );

printf( "The correct result is: %lf\n\n", exp( 1.0 ) );
s[2] = 1.0 - s[0] - s[1];
for( int w = 1; w <= 7; w++ ) {
  printf(
          "%d: %lf   %lf  \n", 
          w, 
          (u.filter( 0, w ))( r ),
          (v.filter( 0, w ))( s ) 
        );
//   cout << w << ":   " << (u.filter( 0, w ))( r ) 
//             << "    " << (v.filter( 0, w ))( s ) 
//             << "\n";
  }
}
End of coordTest.cc
echo evaltest.cc 1>&2
cat >evaltest.cc <<'End of evaltest.cc'
#include <stdlib.h>
#include <stdio.h>
// #include <stream.h>
#include "mxyzptlk.rsc"


//***************************************************************
//***************************************************************
//***************************************************************
// 
//
//	Test program.
//

main() {

printf( "The correct result is: %lf\n\n", exp( 1.0 ) );
// cout << "The correct result is: " << exp( 1.0 ) << "\n\n";

double r[3];
r[0] = 0.5;
r[1] = 0.4;
r[2] = 0.0;
DASetup( 3, 7, r );

DA x, y, z;
DA u, v;

x.setVariable(0);
y.setVariable(1);
z.setVariable(2);

u = exp( x );
v = exp( x + y + z );

double s[3];
r[0] = 1.0;
r[1] = 0.0;
r[2] = 0.0;
s[0] = 0.33;
s[1] = 0.33;
s[2] = 1.0 - s[0] - s[1];

for( int w = 1; w <= 7; w++ ) {
  printf(
          "%d: %lf   %lf  \n", 
          w, 
          (u.filter( 0, w ))( r ),
          (v.filter( 0, w ))( s ) 
        );
//   cout << w << ":   " << (u.filter( 0, w ))( r ) 
//             << "    " << (v.filter( 0, w ))( s ) 
//             << "\n";
  }
}
End of evaltest.cc
echo getTest.cc 1>&2
cat >getTest.cc <<'End of getTest.cc'
#include "mxyzptlk.rsc"


main( int argc, char** argv ) {
static double zero[]   = { 0., 0., 0., 0., 0., 0. }; 
DASetup( 6, 2, 6, zero );                           

DA x, y, z;
DAterm* w;

x.setVariable( 1.0, 0 );
y.setVariable( 2.0, 2 );
x.fixReference();
y.fixReference();

(x*x).peekAt();
(y*y).peekAt();
z = cos( x*x ) + sin( y*y );
z.peekAt();

while ( w = z.get() ) z.peekAt();

}

End of getTest.cc
echo hptest.cc 1>&2
cat >hptest.cc <<'End of hptest.cc'
#include <stdlib.h>
#include <stdio.h>

// ================================================================
// 	Statements needed by any program using DA class


#include "mxyzptlk.rsc"

//***************************************************************
//***************************************************************
//***************************************************************
// 
//
//	Test program.
//

#define	NUMREG 4

#include "menu.hxx"

DA HPreg[ NUMREG ];
DA xda;

int numberOfVar, topWeight;

main() {

ttMenu menu;

void stvr(),
     rvl(),
     rgop(),
     rith(),
     dffr(),
     nlys(),
     clear(),
     eval();

printf("Enter the number of variables and the max order ..> ");
scanf("%d  %d", &numberOfVar, &topWeight );
getchar();
if( numberOfVar > 6 ) {
  printf("Please .. have mercy!!\n");
  printf("I have reduced the number of variables to 6.\n");
  numberOfVar = 6;
  }
if( topWeight > 12 ) {
  printf("Please .. have mercy!!\n");
  printf("I have reduced the maximum order to 12.\n");
  topWeight = 12;
  }
DASetup( numberOfVar, topWeight );

menu.setPrompt( "DA TEST ..> " );
menu.setItem( "Set variable   ", stvr );
menu.setItem( "Peek           ", rvl  );
menu.setItem( "Register op's  ", rgop );
menu.setItem( "Arithmetic     ", rith );
menu.setItem( "Evaluate       ", eval );
menu.setItem( "Differential of", dffr );
menu.setItem( "Analysis       ", nlys );
menu.setItem( "Clear the stack", clear );
menu.go();
}


void stvr(){
double v;
int i;

for( i = NUMREG-1; i > 0; i-- ) HPreg[i] = HPreg[i-1];

printf("Enter the value of the new term ..>  ");
scanf( "%lf", &v );
getchar();
printf("Enter one integer corresponding to the index ..>  ");
scanf( "%d", &i );
getchar();

HPreg[0].setVariable( v, i );

for( i = 0; i < NUMREG; i++ ) HPreg[i].fixReference();
}

void rvl(){
int j = 5;
while( ( j < 0 ) || ( j > NUMREG ) ) {
  printf("Which register?? ..>  ");
  scanf( "%d", &j );
  getchar();
  }
HPreg[j].peekAt();
}

void rgop() {
ttMenu rgopMenu;
void entr(),
     xy(),
     rdown();
rgopMenu.setPrompt( "Reg Ops ..> " );
rgopMenu.setItem( "Set variable   ", stvr  );
rgopMenu.setItem( "Peek           ", rvl   );
rgopMenu.setItem( "Enter          ", entr  );
rgopMenu.setItem( "Interchange    ", xy    );
rgopMenu.setItem( "Rotate         ", rdown );
rgopMenu.go();
}

void rith() {
ttMenu rithMenu;
void neg(),
     plus(),
     minus(),
     mltp(),
     dvd();
rithMenu.setPrompt( "Arithmetic ..> " );
rithMenu.setItem( "Negate         ", neg  );
rithMenu.setItem( "Add            ", plus );
rithMenu.setItem( "Subtract       ", minus );
rithMenu.setItem( "Multiply       ", mltp );
rithMenu.setItem( "Divide         ", dvd  );
rithMenu.go();
}

void eval() {
int i;
double* z;

z = new double[ numberOfVar ];

printf( "Enter %d real numbers ..> ", numberOfVar );
for( i = 0; i < numberOfVar; i++ ) scanf( "%lf ", z+i );
getchar();
printf( "Confirming: You have entered  " );
for( i = 0; i < numberOfVar; i++ ) printf( "%lf ", z[i] );
printf("\n");

getchar();

printf( "\nThe value is %lf.\n", HPreg[0](z)  );

delete z;
}

void dffr() {
int r;
DAterm* p;
printf("Enter the register ..> ");
scanf( "%d", &r );
p = HPreg[r].get();
}

void nlys() {

ttMenu nlysMenu;

void drvt(),
     fltr(),
     xpon(),
     ssnn(),
     ccss(),
     ssqq(),
     llgg();

nlysMenu.setPrompt( "Analysis ..> " );
nlysMenu.setItem( "differentiate  ", drvt );
nlysMenu.setItem( "filter         ", fltr );
nlysMenu.setItem( "exp()          ", xpon );
nlysMenu.setItem( "sin            ", ssnn );
nlysMenu.setItem( "cos            ", ccss );
nlysMenu.setItem( "sqrt           ", ssqq );
nlysMenu.setItem( "log            ", llgg );
nlysMenu.go();

}

void clear() {
int i;
for( i = 0; i < NUMREG; i++ ) HPreg[i].setVariable( 1.0, -3 );
}


void neg(){
DA dummy;
dummy = HPreg[0];
HPreg[0] = - dummy;
}

void plus() {
int i;
HPreg[0] = HPreg[1] + HPreg[0];
for( i = 1; i < NUMREG - 1; i++ ) HPreg[i] = HPreg[i+1];
}

void minus() {
int i;
HPreg[0] = HPreg[1] - HPreg[0];
for( i = 1; i < NUMREG - 1; i++ ) HPreg[i] = HPreg[i+1];
}

void mltp() {
int i;
HPreg[0] = HPreg[1]*HPreg[0];
for( i = 1; i < NUMREG - 1; i++ ) HPreg[i] = HPreg[i+1];
}

void dvd() {
int i;
HPreg[0] = HPreg[1]/HPreg[0];
for( i = 1; i < NUMREG - 1; i++ ) HPreg[i] = HPreg[i+1];
}


void fltr() {
int source, sink, weightLo, weightHi;
printf("Enter source register ..>  ");
scanf( "%d", &source );
getchar();
printf("Enter sink register ..>  ");
scanf( "%d", &sink );
getchar();
printf("Enter desired weight range ..>  ");
scanf( "%d %d", &weightLo, &weightHi );
getchar();
HPreg[ sink ] = HPreg[ source ].filter( weightLo, weightHi );
}

void xpon() {
HPreg[0] = exp( HPreg[0] );
}

void ssnn() {
HPreg[0] = sin( HPreg[0] );
}

void ccss() {
HPreg[0] = cos( HPreg[0] );
}

void ssqq() {
HPreg[0] = sqrt( HPreg[0] );
}

void llgg() {
HPreg[0] = log( HPreg[0] );
}

void drvt() {
int i;
int* n;

n = new int[ numberOfVar ];

printf( "Enter %d integers for the index set ..> ", numberOfVar );
for( i = 0; i < numberOfVar; i++ ) scanf( "%d", n+i );
getchar();
printf( "Confirming: You have entered  " );
for( i = 0; i < numberOfVar; i++ ) printf( "%d ", n[i] );
printf("\n");

HPreg[0] = HPreg[0].D( n );
}


void entr() {
int i;
for( i = NUMREG - 1; i > 0; i-- ) {
  HPreg[i] = HPreg[i-1];
  }
}

void xy() {
DA z;
z = HPreg[1];
HPreg[1] = HPreg[0];
HPreg[0] = z;
}

void rdown() {
DA z;
int i;
z = HPreg[0];
for( i = 0; i < NUMREG - 1; i++ ) HPreg[i] = HPreg[i+1];
HPreg[NUMREG] = z;
}

End of hptest.cc
echo memoryTest.cc 1>&2
cat >memoryTest.cc <<'End of memoryTest.cc'
#include "mxyzptlk.rsc"
#include <stdlib.h>

main( int argc, char** argv )
{
int i, j, count;
DA x, y, z;
static double zero[]   = { 0., 0., 0., 0., 0., 0. }; 
double ff;
DASetup( 2, 2, 2, zero );                           

DA one( 1.0 );
ff = 1.0 / sqrt( 2.0 );
if( argc == 1 ) exit(0);
count = atoi( argv[1] );
i = 0;
x.setVariable( 0.0, 0 );
y.setVariable( 0.0, 2 );
x.fixReference();
y.fixReference();
while( i < count ) {
  for( j = 0; j < 100; j++ ) {
    x.setVariable( 0.0, 0 );
    y.setVariable( 0.0, 2 );
    x.fixReference();
    y.fixReference();
    z = y + x; 
    y = ( y - x );  
//    z = ff*( y + x );       // These two lines are fine separately,
    y = ff*( y - x )*z;     // but disaster together!!!???
//     x = y / ( one + z );
//     x = z;
    i++;
//     printf( "\n=======================================\n" );
  }
  printf( "%d\n", i );
  x.peekAt();
  y.peekAt();
  z.peekAt();
}
}
End of memoryTest.cc
echo pbtest.cc 1>&2
cat >pbtest.cc <<'End of pbtest.cc'
#include <stdlib.h>
#include <stdio.h>
// #include <stream.h>

// ================================================================
// 	Statements needed by any program using DA class

#include "mxyzptlk.rsc"


//***************************************************************
//***************************************************************
//***************************************************************
// 
//
//	Test program.
//

main() {

double u1, u2, v1, v2;
double r[4];
r[0] = u1 = 0.32;
r[1] = u2 = 0.5;
r[2] = v1 = -3.1;
r[3] = v2 = 1.5;

DASetup( 4, 3, 4, r );

double w, y, z, answer;
DA a, b, c;
DA x1, x2, p1, p2;
DA pb;

x1.setVariable(0);
x2.setVariable(1);
p1.setVariable(2);
p2.setVariable(3);

a = (x1*x1) * (x2*x2*x2) * p1 * (p2*p2*p2*p2);
b = sin( x1 * (p2*p2) * (x2*x2*x2) );
pb = a^b;
printf( "Computed by DA: %lf\n", pb.standardPart() );
// cout << "Computed by DA: " << pb.standardPart() << "\n";

w = (u1*u1) * (u2*u2*u2) * v1 * (v2*v2*v2*v2);
y =      u1 * (v2*v2) * (u2*u2*u2)  ;
z = cos( y );
answer = w*y*z*( 6.0/(u2*v2) - 1.0/(u1*v1) - 12.0/(u2*v2) );
printf( "Exact answer  : %lf\n", answer );
// cout << "Exact answer  : " << answer << "\n";

printf( 
      "And also      : %lf\n\n",
      ( ( (x1*x1) * (x2*x2*x2) * p1 * (p2*p2*p2*p2) ) ^
        ( sin( x1 * (p2*p2) * (x2*x2*x2) ) ) 
      ).standardPart()
      );
// cout << "And also      : " 
//      << ( ( (x1*x1) * (x2*x2*x2) * p1 * (p2*p2*p2*p2) ) ^
//           ( sin( x1 * (p2*p2) * (x2*x2*x2) ) ) 
//         ).standardPart()
//      << "\n\n";

// -- Test of the Jacobi identity
c = exp( p1*x1 + p2*x2 );
pb = (a^(b^c)) + (b^(c^a)) + (c^(a^b));
printf( "Jacobi identity\n" );
// cout << "Jacobi identity\n";
pb.peekAt();

}

End of pbtest.cc
echo Lie_A_Test.cc 1>&2
cat >Lie_A_Test.cc <<'End of Lie_A_Test.cc'
/*
**
** Demo program:
** 
** A one-liner: calculates the Lie
** derivative of r**2 under the vector
** field of a damped harmonic oscillator.
** 
** --- Leo Michelotti
** --- August 10, 1993
**
*/

#include <stdlib.h>
#include "mxyzptlk.rsc"

main( int argc, char** argv ) {

 if( argc != 3 ) {
  printf( "\n" );
  printf( "Usage: %s  x  y \n", argv[0] );
  printf( "\n" );
  exit(0);
 }

 DASetup( 2, 1, 2 );

 coord x( atof( argv[1] ) );   
 coord p( atof( argv[2] ) );
 LieOperator v;
 
 v.setVariable(      p, 0 );
 v.setVariable( -x - p, 1 );

 printf( "Answer is: %lf\n", 
         ( v^( x*x + p*p ) ).standardPart() );
}
End of Lie_A_Test.cc
echo Lie_B_Test.cc 1>&2
cat >Lie_B_Test.cc <<'End of Lie_B_Test.cc'
/*
**
** Demo program:
** 
** Uses the exponential map method to 
** integrate the harmonic oscillator
** vector field and obtain the linear
** flow map, or more properly, its
** derivatives.
**
** This particular appliation suffers
** from its being a linear vector field,
** so that the derivatives are not 
** calculated exactly.
** 
** --- Leo Michelotti
** --- August 10, 1993
**
*/

#include "mxyzptlk.rsc"

main() {
 DASetup( 2, 6, 2 );
 coord x( 0.0 );
 coord y( 0.0 );
 LieOperator v;
 v.setVariable(  y, 0 );
 v.setVariable( -x, 1 );
 v.expMap( 6.2831853072, x ) .printCoeffs();
 v.expMap( 6.2831853072, y ) .printCoeffs();
}
End of Lie_B_Test.cc
echo Lie_C_Test.cc 1>&2
cat >Lie_C_Test.cc <<'End of Lie_C_Test.cc'
/*
**
** Demo program:
** 
** Uses the exponential map method to 
** integrate x' = x.
** 
** --- Leo Michelotti
** --- August 10, 1993
**
*/

#include <stdlib.h>
#include "mxyzptlk.rsc"

main( int argc, char** argv ) {

 if( argc != 3 ) {
  printf( "\n" );
  printf( "Usage: %s  d  t                           \n", argv[0] );
  printf( "where  d  = degree of highest derivative  \n" );
  printf( "       t  = upper limit of integration.   \n" );
  printf( "\n" );
  exit(0);
 }

 int    deg   = atoi( argv[1]  );
 double limit = atof( argv[2]  );

 DASetup( 2, deg, 1 );
 coord x( 0.0 ), t( limit );
 LieOperator v;
 v.setVariable( x, 0 );

 printf( "\n-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-\n" );
 printf( "Upper limit is double variable.\n" );
 v.expMap( limit, x ) .printCoeffs();
 printf( "\n-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-\n" );
 printf( "Upper limit is DA variable.\n" );
 v.expMap( t, x ) .printCoeffs();
 
}
End of Lie_C_Test.cc
echo Lie_D_Test.cc 1>&2
cat >Lie_D_Test.cc <<'End of Lie_D_Test.cc'
/*
**
** Demo program:
** 
** Uses the exponential map method to 
** integrate a non-Hamltonian, 
** nonlinear vector
** field that spirals inward.
** 
** --- Leo Michelotti
** --- August 10, 1993
**
*/

#include <stdlib.h>
#include "mxyzptlk.rsc"

main( int argc, char** argv ) {

 if( argc != 3 ) {
  printf( "\n" );
  printf( "Usage: %s  d  t                           \n", argv[0] );
  printf( "where  d  = highest degree derivative     \n" );
  printf( "       t  = upper limit of integration.   \n" );
  printf( "\n" );
  exit(0);
 }

 int    deg   = atoi( argv[1] );
 double limit = atof( argv[2] );

 DASetup( 3, deg, 2 );
 coord x( 0.0 );
 coord y( 0.0 );
 coord t( limit );
 LieOperator v;
 DA f, g;
 
 f = 1.0 - exp( - x*x - y*y );
 v.setVariable(          f*y, 0 );
 v.setVariable( -f*( x + y ), 1 );

 printf( "\n-*-*-*-  x  coefficients  -*-*-*-*-*\n" );
 ( f = v.expMap( t, x ) ).printCoeffs();
 printf( "\n-*-*-*-  y  coefficients  -*-*-*-*-*\n" );
 ( g = v.expMap( t, y ) ).printCoeffs();

 // Enter a point and evaluate the map.
 double s[3];
 double tlo, thi, delta;
 printf( "Enter x and y:  ");
 scanf( "%lf  %lf", s, s+1);
 getchar();
 printf( "Enter range and increment  ");
 scanf( "%lf  %lf  %lf", &tlo, &thi, &delta );
 getchar();
 s[2] = limit;
 printf( "( %lf, %lf ) maps to ( %lf, %lf )\n\n",
          s[0], s[1] , f(s), g(s) );
 for( limit = tlo; limit < thi + delta / 2.0 ; limit += delta ) {
  s[2] = limit;
  printf( "%lf: ( %lf, %lf )\n",
          limit, f(s), g(s) );
 }
}
End of Lie_D_Test.cc
echo Lie_E_Test.cc 1>&2
cat >Lie_E_Test.cc <<'End of Lie_E_Test.cc'
/*
**
** Demo program:
** 
** Uses the exponential map method to 
** integrate equations of motion of a 
** sextupolar type vector field.
** 
** As part of demo, we introduce the
** class field to make writing the
** vector field easier, as well as to
** illustrate style and elegance.
** 
** --- Leo Michelotti
** --- August 10, 1993
**
*/

#include <stdlib.h>
#include "mxyzptlk.rsc"

class field {
 private:
  DA v[3];
 public:
  field( coord&, coord&, coord& );
  DA operator[]( int );
};

field::field( coord& x, coord& y, coord& z ) {
 v[0] = - 2.0*x*y;
 v[1] =   x*x - y*y;
 v[2] =   0.0;
}

DA field::operator[]( int i ) {
 return v[i];
}

main( int argc, char** argv ) {

 if( argc != 4 ) {
  printf( "\n" );
  printf( "Usage: %s  d  m  t                        \n", argv[0] );
  printf( "where  d  = highest order derivative      \n" );
  printf( "       m  = mass                          \n" );
  printf( "       t  = upper limit of integration.   \n" );
  printf( "\n" );
  exit(0);
 }

 int    deg = atoi( argv[1] );
 double msq = atof( argv[2] );
 double t   = atof( argv[3] );

 msq = msq*msq;
 DASetup( 6, deg, 6 );

 coord  x( 0.0 ),  y( 0.0 ),  z( 0.0 );
 coord px( 0.0 ), py( 0.0 ), pz( 0.0 );
 DA E = sqrt( px*px + py*py + pz*pz + msq );
 field B( x, y, z );

 LieOperator v;
 v.setVariable(  px		       / E, 0 );
 v.setVariable(  py		       / E, 1 );
 v.setVariable(  pz		       / E, 2 );
 v.setVariable(  ( py*B[2] - pz*B[1] ) / E, 3 );
 v.setVariable(  ( pz*B[0] - px*B[2] ) / E, 4 );
 v.setVariable(  ( px*B[1] - py*B[0] ) / E, 5 );

 printf( "\n-*-*-*-  x  coefficients  -*-*-*-*-*\n" );
 v.expMap( t, x ).printCoeffs();
 printf( "\n-*-*-*-  y  coefficients  -*-*-*-*-*\n" );
 v.expMap( t, y ).printCoeffs();
 printf( "\n-*-*-*-  z  coefficients  -*-*-*-*-*\n" );
 v.expMap( t, z ).printCoeffs();
}
End of Lie_E_Test.cc
echo Limit_B_Test.cc 1>&2
cat >Limit_B_Test.cc <<'End of Limit_B_Test.cc'
/*
**
** Demo program:
** Calculates the arctangent function and 
** derivatives at a point specified on the
** command line.
** --- Leo Michelotti
** --- February 12, 1993
** 
*/

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "mxyzptlk.rsc"

main( int argc, char** argv ) {

 const int  dim       = 1;

 if( argc != 4 ) {
  printf( "Usage: Limit_B_Test <point of evaluation> <iterations> <max order>\n" );
  exit(0);
 }

 double  x          = atof( argv[1] );
 int     iterations = atoi( argv[2] );
 int     maxWeight  = atoi( argv[3] );
 DASetup( dim, maxWeight, dim );

 double answer = atan(x);
 printf ( "atan( %lf ) = %le\n", x, answer );

 coord m( x );

 DA u, s, c;
 int d[1];
 d[0] = 0;
 u = m;

 for( int i = 0; i < iterations; i++ ){
  s = sin( u );
  c = cos( u );
  u = u - c*( s - m*c );
  printf( "%le\n", u.standardPart() );
  for( int j = 1; j < 6; j++ ) {
   d[0] = j;
   printf( "%le  ", u.derivative(d) );
  }
  printf( "\n" );
 }

}
End of Limit_B_Test.cc
echo Limit_C_Test.cc 1>&2
cat >Limit_C_Test.cc <<'End of Limit_C_Test.cc'
/*
**
** Demo program:
** Calculates the DA expansion of 
** arctan( sqrt( y^2 + z^2 )) at
** a point of evaluation specified 
** command line.  Ten iterations of
** a Newton's method algorithm are
** used.
**
** --- Leo Michelotti
** --- February 12, 1993
** 
*/

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include "mxyzptlk.rsc"

main( int argc, char** argv ) {

 if( argc != 3 ) {
  printf( "Usage: Limit_C_Test <y> <z> <degree>\n" );
  exit(0);
 }

 double ay        = atof( argv[1] );
 double az        = atof( argv[2] );
 int    degree    = atoi( argv[3] );

 DASetup( 2, degree, 2 );

 coord y( ay );
 coord z( az );

 DA m, u, s, c;
 m = sqrt( y*y + z*z );
 u = m;

 for( int i = 0; i < 10; i++ ){
  s = sin( u );
  c = cos( u );
  u = u - c*( s - m*c );
  printf( "%le\n", u.standardPart() );
  u.peekAt();
  printf( "\n" );
 }

}

/*
The output was tested against Mathematica. 
Here is the result:

churchy 2: math

Mathematica 2.0 for SPARC
Copyright 1988-91 Wolfram Research, Inc.

In[1]:= Series[ArcTan[Sqrt[(1+y)*(1+y) + z*z]], {y,0,5}, {z,0,5}]

              2    4                 2      4
        Pi   z    z        6    1   z    7 z        6
Out[1]= -- + -- - -- + O[z]  + (- - -- + ---- + O[z] ) y + 
        4    4    8             2   2     16
 
                2       4                        2       4
        1    5 z    29 z        6   2    1    5 z    47 z        6   3
>    (-(-) + ---- - ----- + O[z] ) y  + (-- - ---- + ----- + O[z] ) y  + 
        4     8      32                  12    8      32
 
         2       4                         2        4
      9 z    67 z        6   4      1     z    179 z        6   5       6
>    (---- - ----- + O[z] ) y  + (-(--) - -- + ------ + O[z] ) y  + O[y]
       16     32                    40    2      64

In[2]:= ExpandAll[%]

              2    4                 2      4
        Pi   z    z        6    1   z    7 z        6
Out[2]= -- + -- - -- + O[z]  + (- - -- + ---- + O[z] ) y + 
        4    4    8             2   2     16
 
                2       4                        2       4
        1    5 z    29 z        6   2    1    5 z    47 z        6   3
>    (-(-) + ---- - ----- + O[z] ) y  + (-- - ---- + ----- + O[z] ) y  + 
        4     8      32                  12    8      32
 
         2       4                         2        4
      9 z    67 z        6   4      1     z    179 z        6   5       6
>    (---- - ----- + O[z] ) y  + (-(--) - -- + ------ + O[z] ) y  + O[y]
       16     32                    40    2      64

In[3]:= Quit[];

*/
End of Limit_C_Test.cc
echo Limit_D_Test.cc 1>&2
cat >Limit_D_Test.cc <<'End of Limit_D_Test.cc'
/*
**
** Demo program:
** Calculates the DA expansion of arcsin 
** using Newton's method.
** --- Leo Michelotti
** --- February 12, 1993
** 
*/

#include <stdlib.h>
#include <stdio.h>
#include "mxyzptlk.rsc"

main( int argc, char** argv ) {

 if( argc != 3 ) {
  printf( "Usage: Limit_D_Test <degree> <iterations>\n" );
  exit(0);
 }

 int deg = atoi( argv[1] );
 int n   = atoi( argv[2] );
 DASetup( 1, deg );

 coord x(0.0);
 DA u;

 u = x;
 for( int i = 0; i < n; i++ ) 
  u -= ( sin(u) - x ) / cos(u);
 u.printCoeffs();
}
End of Limit_D_Test.cc
echo Limit_E_Test.cc 1>&2
cat >Limit_E_Test.cc <<'End of Limit_E_Test.cc'
/*
**
** Demo program:
** Calculates the DA expansion of the 
** inverse of an "arbitrary" function
** - in this case, x + x^2 + x^3 - 
** using Newton's method.
** --- Leo Michelotti
** --- February 12, 1993
** 
*/

#include <stdlib.h>
#include <stdio.h>
#include "mxyzptlk.rsc"

main( int argc, char** argv ) {

 if( argc != 3 ) {
  printf( "Usage: Limit_E_Test <degree> <iterations>\n" );
  exit(0);
 }

 int deg = atoi( argv[1] );
 int n   = atoi( argv[2] );
 DASetup( 1, deg, 1 );

 coord x(0.0);
 DA u;
 DA f, g;
 int d = 1;

 f = x + x*x + x*x*x;
 g = f.D(&d);
 
 u = x;
 for( int i = 0; i < n; i++ ) {
  u -= ( f(&u) - x ) / g(&u);
  u.peekAt();
 }

 f(&u).peekAt();
}
End of Limit_E_Test.cc
echo Inverse_A_Test.cc 1>&2
cat >Inverse_A_Test.cc <<'End of Inverse_A_Test.cc'
/*
**
** Test of DA::inverse() function.
** 
** --- Leo Michelotti
** --- September 16, 1993
**
*/

#include <stdlib.h>
#include <stdio.h>
#include "mxyzptlk.rsc"

main( int argc, char** argv ) {

 if( argc != 2 ) {
  printf( "Usage: %s <degree>\n", argv[0] );
  exit(0);
 }

 int deg = atoi( argv[1] );
 DASetup( 1, deg, 1 );

 coord x(0.0);
 DA f, g;

 f = x + x*x + x*x*x;
 g = f.inverse();
 g.printCoeffs();
 f(&g).printCoeffs();
 g(&f).printCoeffs();
}
End of Inverse_A_Test.cc
echo Inverse_B_Test.cc 1>&2
cat >Inverse_B_Test.cc <<'End of Inverse_B_Test.cc'
/*
**
** Test of LieOperator inversion algorithm.
** 
** --- Leo Michelotti
** --- September 15, 1993
**
*/

#include <stdlib.h>
#include <stdio.h>
#include "mxyzptlk.rsc"

main( int argc, char** argv ) {
 
 if( argc != 2 ) {
  printf( "\n" );
  printf( "Usage: %s  limit                          \n", argv[0] );
  printf( "\n" );
  exit(0);
 }

 int limit = atoi( argv[1] );

 DASetup( 3, 4, 3 );

 coord        x(0.0), y(0.0), z(0.0);
 LieOperator  u, v, w;
 LieOperator  id( "ident" );

 w.setVariable( x + 3.0*y + x*y - y*z,   0 );
 w.setVariable( y - x + z + x*z + y*y,   1 );
 w.setVariable( z + 2.0*x + y*z - x*y*z, 2 );

 doublematrix M( 3, 3 );
 M = w.Jacobian();
 M = M.inverse();

 u = M*id;
 v = ( w(u) - id );
 int i = 0;
 while ( ( i++ < limit ) && ( v != 0.0 ) ) {
  u = u - M*v;
  v = ( w(u) - id );
 }

 printf( "\n====== u.printCoeffs(); =================\n" );
 u.printCoeffs();

 printf( "\n====== w.printCoeffs(); =================\n" );
 w.printCoeffs();

 printf( "\n====== w(u).printCoeffs(); =================\n" );
 w(u).printCoeffs();

 printf( "\n====== u(w).printCoeffs(); =================\n" );
 u(w).printCoeffs();

 if( i >= limit ) printf( "++++ Limit exceeded! ++++\n" );
}
End of Inverse_B_Test.cc
echo Inverse_C_Test.cc 1>&2
cat >Inverse_C_Test.cc <<'End of Inverse_C_Test.cc'
/*
**
** Test of LieOperator::inverse() 
** This does not work well because of retaining
** negligible terms (of order 10**-16).  I do
** not yet understand why these terms are being
** kept.
** 
** --- Leo Michelotti
** --- September 15, 1993
**
*/

#include <stdlib.h>
#include <stdio.h>
#include "mxyzptlk.rsc"

main( int argc, char** argv ) {
 
 DASetup( 3, 4, 3 );

 coord        x(0.0), y(0.0), z(0.0);
 LieOperator  w;
 LieOperator  u;

 w.setVariable( x + 3.0*y + x*y - y*z,   0 );
 w.setVariable( y - x + z + x*z + y*y,   1 );
 w.setVariable( z + 2.0*x + y*z - x*y*z, 2 );

 u = w.inverse();

 printf( "\n====== w.printCoeffs(); =================\n" );
 w.printCoeffs();

 printf( "\n====== u.printCoeffs(); =================\n" );
 u.printCoeffs();

 printf( "\n====== w(u).printCoeffs(); =================\n" );
 w(u).printCoeffs();

 printf( "\n====== u(w).printCoeffs(); =================\n" );
 u(w).printCoeffs();
}
End of Inverse_C_Test.cc
echo newtonAA.cc 1>&2
cat >newtonAA.cc <<'End of newtonAA.cc'
/* 
 * June 14, 1991
 * 
 * Serves as a model for a general purpose tool
 * for finding the DA representation of implicitly 
 * defined functions.  In this case we want to 
 * find x(m) such that 
 * 
 *   x(m) = cos( m*x(m) )
 * 
 * -- Leo Michelotti
 * 
 */

#include <stdio.h>
#include "mxyzptlk.rsc"

// ...... Globals ..........................................
const  double tolerance = 1.0e-8;
static int    dfNdx[]   = { 0, 1 };
const  int    maxIter   = 15;

double  m;
DA      G, md, xd;

// ...... Functions ........................................
DA     F( DA& x )    { return ( x - cos( md*x ) ); }
double F( double x ) { return ( x - cos( m *x ) ); }

double solve( double m, double x ) {
double xOld, s[2];
int    i;

md.setVariable( m, 0 );
xd.setVariable( x, 1 );
md.fixReference();
xd.fixReference();

s[0] = m;    i = 0;
s[1] = x;    xOld = 123456789.0;

while( ( tolerance < fabs( x - xOld ) ) 
   &&  ( i++ < maxIter ) 
     ) {
  xOld = x;
  G = F( xd );
  G = xd - ( G / G.D( dfNdx ) ); // Sets up Newton's algorithm.
  x = G( s );                    // Recursive step.
  md.setVariable( m, 0 );
  xd.setVariable( x, 1 );
  md.fixReference();
  xd.fixReference();
  s[1] = x;
}
G = F( xd );                     // This is done to assure that G has
G = xd - ( G / G.D( dfNdx ) );   // a reference point identical to
                                 // md and xd, at least at the start.
return x;
}

// ...... Main program .....................................
main() {
const int  dim       = 2;
const int  maxWeight = 5;
double     x;
double     mLo, mHi, step, s[2];
int        i;

DASetup( dim, maxWeight, dim ); // This form required
DA sd[2], Y;                    // to do multiEval( DA* )

printf("Enter the value of m\n");
printf("and an initial guess for x: ");
scanf("%lf %lf", &m, &x );
getchar();
printf( "You have chosen m = %lf and x = %lf\n", m, x );

// ...... First solve for the value of x(m) ................
x = solve( m, x );
printf( "x( %lf ) = %lf   ", m, x );
printf( "Defect is: %lf\n", F(x) );

// ...... Now calculate the derivatives ....................
double krunge[2];
sd[0] = md;
sd[1] = xd;
for( i = 0; i < maxIter; i++ ) {
  xd = G( sd );
  sd[0] = md;
  sd[1] = xd;
  for( int j = 0; j < 2; j++ ) 
    krunge[j] = sd[j].standardPart();   // In a perfect world this
  G.fixReference( krunge );             // would not be necessary.
}

xd.peekAt();

// ...... And validate .....................................
printf( "\nThe real test:\n" );
printf( "Enter mLo, mHi, and step: " );
scanf("%lf %lf %lf", &mLo, &mHi, &step );
getchar();

Y = xd;          // xd is a global variable, and solve() will 
                 // change its value.
s[1] = 1.0;      // The value of s[1] no longer matters.
for( m = mLo; m < (mHi + 0.5*step); m+=step ){
  x = solve( m, x );
  s[0] = m;  
  printf( "%lf  %lf  %lf  %lf\n",
          m,   x,   Y( s ),   F(x)
        );
}

} // End function: main()
End of newtonAA.cc
