using math, reduce;

// Initialization (this requires the REDUCE image file in the current dir).
reduce::start "/usr/local/pure-lang/pure-reduce/reduce.img" {"-v"};

// Verbosity level. 0 means no noise at all.
reduce::verbosity 0;


// Pauli matrices (sigma 1..3).
// see e.g. http://en.wikipedia.org/wiki/Pauli_matrices
let s0 = {1,0;0,1} ; 
let s1 = {0,1;1,0} ;
let s2 = {0,-i;i,0};
let s3 = {1,0;0,-1};

// Check the identities: s1^2=s2^2=s3^2=-i*s1*s2*s3=s0
// Note: s1^2 or s1*s1 works here.
let r1 = simplify $ (s1*s1) ;
let r2 = simplify $ (s2*s2) ;
let r3 = simplify $ (s3*s3) ;
let r4 = simplify $ (-i*s1*s2*s3) ;
let r5 = all (\s->s==s0) [r1,r2,r3,r4] ; // -> 1

// All det s_i = -1, i=1,2,3 (obvioulsy:)
let r6 = map (simplify . det) [s1,s2,s3] ; // -> [-1,-1,-1]

// Eigenvalues/-vectors of s2
let r7 = simplify $ mateigen s2 q; 
let r8 = map head r7; // -> [q-1,q+1] => Eigenvalues q=+/-1
let r9 = map (head .tail) r7 ; // [1,1], multiplicities
let r10 = map last r7 ; // [c1*{-i,1},c2*{i,1}], eigenvectors

// Transpose
let r11 = map (simplify .tp) [s1,s2,s3] ; // -> [s1',s2',s3']  

// Trace 
let r11 = map (simplify .trace) [s1,s2,s3] ; // -> [tr s1,tr s2,tr s3]=[0,0,0]

// Cofactor 
let r12 = simplify $ cofactor s2 1 1 ; // -> 0.

// Nullspace of s2+{0,i;0,0}
let r13 = simplify $ (s2+{0,i;0,0}) ; // -> {0,0;i,0}

// Rank
let r14 = map (simplify . rank) [s0,s1,s2,s3] ; // -> [2,2,2,2] 

// Inverse  (simply 1/matrix)
let r15 = simplify $ 1/s2 ; // -> {0,1/i;(-1)/i,0}
let r16 = simplify $ s2*r15 ; // -> {1,0;0,1} == s0 as it should be.

// Solving without "solve":
// a11*x(1) + a12*x(2) = y1
// a21*x(1) + a22*x(2) = y2
#!
let r17 = simplify $ (1/{a11,a12;a21,a22}*{y1;y2}) ; // A^-1 * y' ;
#!
// -> {(-a12*y2+a22*y1)/(a11*a22-a12*a21);(a11*y2-a21*y1)/(a11*a22-a12*a21)}

////////////
// ALGINT //
////////////
let r18 = simplify $ intg  (sqrt(x+sqrt(x^2+1))/x) x ; 
// -> intg (sqrt(sqrt (x^2+1)+x)/x) x, i.e. unevaluated !
// load the ALGINT package and retry:
reduce::load "algint" ;
let r19 = simplify $ intg  (sqrt(x+sqrt(x^2+1))/x) x ; 
// -> atan ((sqrt (sqrt (x^2+1)+x)*sqrt (x^2+1)-sqrt (sqrt (x^2+1)+x)*x-sqrt 
// (sqrt (x^2+1)+x))/2)+2*sqrt (sqrt (x^2+1)+x)+log 
// (sqrt (sqrt (x^2+1)+x)-1)-log (sqrt (sqrt (x^2+1)+x)+1)

////////////
// LIMITS //
////////////
let r20 = simplify $ limit (x*sin(1/x)) x (quote infinity) ; // -> 1
let r21 = simplify $ limit (1/x) x 0 ; // -> infinity
// Notes:
// - apparently it's not necessary to do a package "load",
// - map (Pure)inf <-> (Reduce)infinity ?


//////////////
// ODESOLVE //
//////////////
let r22 = lisp ('depend [y,x]); // declare: y depends on x
let r23 = simplify $ odesolve [df y x == x^2+(quote e)^x] [y] x;
let r24 = simplify $ odesolve [(df y x 2) == y] [y] x [[x==0,y==A],[x==1,y==B]] ;
//
// r23: y' = x^2 + e^x
// r24: y'' = y with boundary cond: y(0)=A, y(1)=B
// see manual for initial cond and other features
//


//////////////////////////////////////////
// SUM: A package for series summation ///
//////////////////////////////////////////
// This package implements the Gosper algorithm for the summation of series.
let r25 = 0;//simplify $ sum n^3 n 1 N ;
let r26 = 0;//simplify $ sum (a+k*r) k 0 (n-1) ;
let r27 = 0;//simplify $ sum (1/((p+(k-1)*q)*(p+k*q))) k 1 (n+1) ;
let r28 = 0;//simplify $ prod (k/(k-2)) k 1 N ;
//***** Declare min=>min operator ?  (Y or N)

//////////////////////////////////////////
//TAYLOR: Manipulation of Taylor series //
//////////////////////////////////////////
def exp x = (quote e)^x;
let r29 = simplify $ taylor (exp(x^2+y^2)) x 0 2 y 0 2 ;
//** taylor_ [[[0],[0]]:1:1,[[0],[2]]:1:1,[[2],[0]]:1:1,[[2],[2]]:1:1]
//           [[[x],0,2,3],[[y],0,2,3]] () ()
let r30 = simplify $ taylor ((quote e)^x) x 0 3;
//** taylor_ [[[0]]:1:1,[[1]]:1:1,[[2]]:1:2,[[3]]:1:6] [[[x],0,3,4]] () ()

// should give:
/*
Reduce (Free CSL version), 14-Apr-11 ...

1: taylor(e^(x^2+y^2),x,0,2,y,0,2);

     2    2    2  2      3  3
1 + y  + x  + y *x  + O(x ,y )

2: taylor(e^x,x,0,3);

         1   2    1   3      4
1 + x + ---*x  + ---*x  + O(x )
         2        6

3:
*/