yacas 1.3.6-2+b1 / usr / share / yacas / scripts / sums.rep / code.ys

/*  this is disabled because some functions seem to implicitly define Min / Max with a different number of args, and then Yacas is confused if it hasn't loaded all the Function() declarations beforehand.
FIXME
/// Min, Max with many arguments
*/
Retract("Min", 1);
Retract("Min", 2);
Retract("Min", 3);
Retract("Max", 1);
Retract("Max", 2);
Retract("Max", 3);
//Function() Min(list);
//Function() Max(list);

//Function() Min(l1, l2);
//Function() Max(l1, l2);
Function() Min(l1, l2, l3, ...);
Function() Max(l1, l2, l3, ...);

10 # Min(_l1, _l2, l3_IsList) <-- Min(Concat({l1, l2}, l3));
20 # Min(_l1, _l2, _l3) <-- Min({l1, l2, l3});

10 # Max(_l1, _l2, l3_IsList) <-- Max(Concat({l1, l2}, l3));
20 # Max(_l1, _l2, _l3) <-- Max({l1, l2, l3});
/**/
10 # Min(l1_IsList,l2_IsList) <-- Map("Min",{l1,l2});
10 # Max(l1_IsList,l2_IsList) <-- Map("Max",{l1,l2});

20 # Min(l1_IsRationalOrNumber,l2_IsRationalOrNumber) <-- If(l1<l2,l1,l2);
20 # Max(l1_IsRationalOrNumber,l2_IsRationalOrNumber) <-- If(l1>l2,l1,l2);

30 # Min(l1_IsConstant,l2_IsConstant) <-- If(N(Eval(l1-l2))<0,l1,l2);
30 # Max(l1_IsConstant,l2_IsConstant) <-- If(N(Eval(l1-l2))>0,l1,l2);

// Min and Max on empty lists
10 # Min({}) <-- Undefined;
10 # Max({}) <-- Undefined;

20 # Min(list_IsList) <--
[
  Local(result);
  result:= list[1];
  ForEach(item,Tail(list)) result:=Min(result,item);
  result;
];
20 # Max(list_IsList) <--
[
  Local(result);
  result:= list[1];
  ForEach(item,Tail(list)) result:=Max(result,item);
  result;
];

30 # Min(_x) <-- x;
30 # Max(_x) <-- x;

/* Factorials */

10 # 0! <-- 1;
10 # (Infinity)! <-- Infinity;
20 # ((n_IsPositiveInteger)!) <-- [
	Check(n <= 65535, "Factorial: Error: the argument " : ( ToString() Write(n) ) : " is too large, you may want to avoid exact calculation");
	MathFac(n);
];

25 # ((x_IsConstant)!)_(FloatIsInt(x) And x>0) <-- (Round(x)!);

30 # ((x_IsNumber)!)_InNumericMode() <-- Internal'GammaNum(x+1);

40 # (n_IsList)! <-- MapSingle("!",n);

/* formulae for half-integer factorials:

(+(2*z+1)/2)! = Sqrt(Pi)*(2*z+1)! / (2^(2*z+1)*z!) for z >= 0
(-(2*z+1)/2)! = Sqrt(Pi)*(-1)^z*z!*2^(2*z) / (2*z)! for z >= 0

Double factorials are more efficient:
	(2*n-1)!! := 1*3*...*(2*n-1) = (2*n)! / (2^n*n!)
	(2*n)!! := 2*4*...*(2*n) = 2^n*n!

*/
/* // old version - not using double factorials
HalfIntegerFactorial(n_IsOdd) _ (n>0) <--
	Sqrt(Pi) * ( n! / ( 2^n*((n-1)/2)! ) );
HalfIntegerFactorial(n_IsOdd) _ (n<0)  <--
	Sqrt(Pi) * ( (-1)^((-n-1)/2)*2^(-n-1)*((-n-1)/2)! / (-n-1)! );
*/
// new version using double factorials
HalfIntegerFactorial(n_IsOdd) _ (n>0) <--
	Sqrt(Pi) * ( n!! / 2^((n+1)/2) );
HalfIntegerFactorial(n_IsOdd) _ (n<0)  <--
	Sqrt(Pi) * ( (-1)^((-n-1)/2)*2^((-n-1)/2) / (-n-2)!! );
//HalfIntegerFactorial(n_IsOdd) _ (n= -1)  <-- Sqrt(Pi);

/* Want to also compute (2.5)! */
40 # (n_IsRationalOrNumber)! _(Denom(Rationalize(n))=2) <-- HalfIntegerFactorial(Numer(Rationalize(n)));

/// partial factorial
n1_IsRationalOrNumber *** n2_IsRationalOrNumber <--
[
	Check(n2-n1 <= 65535, "Partial factorial: Error: the range " : ( ToString() Write(n2-n1) ) : " is too large, you may want to avoid exact calculation");
	If(n2-n1<0,
		1,
		Factorial'partial(n1, n2)
	);
];

/// recursive routine to evaluate "partial factorial" a*(a+1)*...*b
// TODO lets document why the >>1 as used here is allowed (rounding down? What is the idea behind this algorithm?)
2# Factorial'partial(_a, _b) _ (b-a>=4) <-- Factorial'partial(a, a+((b-a)>>1)) * Factorial'partial(a+((b-a)>>1)+1, b);
3# Factorial'partial(_a, _b) _ (b-a>=3) <-- a*(a+1)*(a+2)*(a+3);
4# Factorial'partial(_a, _b) _ (b-a>=2) <-- a*(a+1)*(a+2);
5# Factorial'partial(_a, _b) _ (b-a>=1) <-- a*(a+1);
6# Factorial'partial(_a, _b) _ (b-a>=0) <-- a;


/* Binomials -- now using partial factorial for speed */
// Bin(n,m) = Bin(n, n-m)
10 # Bin(0,0) 		<-- 1;
10 # Bin(n_IsPositiveInteger,m_IsNonNegativeInteger)_(2*m <= n) <-- ((n-m+1) *** n) / m!;
15 # Bin(n_IsPositiveInteger,m_IsNonNegativeInteger)_(2*m > n And m <= n) <-- Bin(n, n-m);
20 # Bin(n_IsInteger,m_IsInteger) <-- 0;

/// even/odd double factorial: product of even or odd integers up to n
1# (n_IsPositiveInteger)!! _ (n<=3) <-- n;
2# (n_IsPositiveInteger)!! <--
[
	Check(n<=65535, "Double factorial: Error: the argument " : ( ToString() Write(n) ) : " is too large, you may want to avoid exact calculation");
	Factorial'double(2+Mod(n, 2), n);
];
// special cases
3# (_n)!! _ (n= -1 Or n=0)<-- 1;

// the purpose of this mess "Div(a+b,2)+1+Mod(Div(a+b,2)+1-a, 2)" is to obtain the smallest integer which is >= Div(a+b,2)+1 and is also odd or even when a is odd or even; we need to add at most 1 to (Div(a+b,2)+1)
2# Factorial'double(_a, _b) _ (b-a>=6) <-- Factorial'double(a, Div(a+b,2)) * Factorial'double(Div(a+b,2)+1+Mod(Div(a+b,2)+1-a, 2), b);
3# Factorial'double(_a, _b) _ (b-a>=4) <-- a*(a+2)*(a+4);
4# Factorial'double(_a, _b) _ (b-a>=2) <-- a*(a+2);
5# Factorial'double(_a, _b) <-- a;

/// double factorial for lists is threaded
30 # (n_IsList)!! <-- MapSingle("!!",n);



/* Sums */

RuleBase("Sum",{sumvar'arg,sumfrom'arg,sumto'arg,sumbody'arg}); 

5  # Sum(_sumvar,sumfrom_IsNumber,sumto_IsNumber,_sumbody)_(sumfrom>sumto) <-- 0;

10 # Sum(_sumvar,sumfrom_IsNumber,sumto_IsNumber,_sumbody)_(sumto<sumfrom) <--
     ApplyPure("Sum",{sumvar,sumto,sumfrom,sumbody});
20 # Sum(_sumvar,sumfrom_IsNumber,sumto_IsNumber,_sumbody) <--
LocalSymbols(sumi,sumsum)[
   Local(sumi,sumsum);
   sumsum:=0;
   For(sumi:=sumfrom,sumi<=sumto,sumi++)
       [
        MacroLocal(sumvar);
        MacroSet(sumvar,sumi);
        sumsum:=sumsum+Eval(sumbody);
       ];
   sumsum;
];

UnFence("Sum",4);
HoldArg("Sum",sumvar'arg);
HoldArg("Sum",sumbody'arg);

Function() Add(val, ...);

10 # Add({}) <-- 0;
20 # Add(values_IsList) <--
[
   Local(i, sum);
   sum:=0;
   ForEach(i, values) [ sum := sum + i; ];
   sum;
];

// Add(1) should return 1
30 # Add(_value) <-- value;


Function("Factorize",{sumvar,sumfrom,sumto,sumbody})
[
   Local(sumi,sumsum);
   sumsum:=1;
   For(sumi:=sumfrom,sumi<=sumto And sumsum!=0,sumi++)
       [
        MacroLocal(sumvar);
        MacroSet(sumvar,sumi);
        sumsum:=sumsum*Eval(sumbody);
       ];
   sumsum;
];
UnFence("Factorize",4);
HoldArg("Factorize",sumvar);
HoldArg("Factorize",sumbody);

Factorize(sumlist_IsList) <--
[
   Local(sumi,sumsum);
   sumsum:=1;
   ForEach(sumi,sumlist)
   [
     sumsum:=sumsum*sumi;
   ];
   sumsum;
];


/*COMMENT FROM AYAL: Jitse, I added some code to make Taylor2 work in the most general case too I believe.
  Could you check to see if you agree with my changes? If that is correct, perhaps we can start calling Taylor2
  by default in stead of Taylor1.
 */
Function("Taylor",{taylorvariable,taylorat,taylororder,taylorfunction})
  Taylor1(taylorvariable,taylorat,taylororder)(taylorfunction);

/*COMMENT FROM AYAL: this is the old slow but working version of Taylor series expansion. Jitse wrote a
 * faster version which resides in taylor.ys, and uses lazy power series. This slow but correct version is still
 * useful for tests (the old and the new routine should yield identical results).
 */
Function("Taylor1",{taylorvariable,taylorat,taylororder,taylorfunction})
[
  Local(n,result,dif,polf);
  [
    MacroLocal(taylorvariable);
    [
      MacroLocal(taylorvariable);
      MacroSet(taylorvariable, taylorat);
      result:=Eval(taylorfunction);
    ];
    If(result=Undefined,
    [
      result:=Apply("Limit",{taylorvariable,taylorat,taylorfunction});
    ]);
/*
    MacroSet(taylorvariable,taylorat);
    result:=Eval(taylorfunction);
*/
  ];
  dif:=taylorfunction;
  polf:=(taylorvariable-taylorat);
  For(n:=1,result != Undefined And n<=taylororder,n++)
  [
    dif:= Deriv(taylorvariable) dif;
    Local(term);
    MacroLocal(taylorvariable);
    [
      MacroLocal(taylorvariable);
      MacroSet(taylorvariable, taylorat);
      term:=Eval(dif);
    ];
    If(term=Undefined,
    [
      term:=Apply("Limit",{taylorvariable,taylorat,dif});
    ]);
      
    result:=result+(term/(n!))*(polf^n);
/*    result:=result+Apply("Limit",{taylorvariable,taylorat,(dif/(n!))})*(polf^n); */
/*
    MacroSet(taylorvariable,taylorat);
    result:=result+(Eval(dif)/(n!))*(polf^n);
*/    
  ];
  result;
];


Function("Subfactorial",{n}) 
[
	n! * Sum(k,0,n,(-1)^(k)/k!);
];

30 # Subfactorial(n_IsList) <-- MapSingle("Subfactorial",n);

// Attempt to Sum series 

LocalSymbols(c,d,expr,from,to,summand,sum,predicate,k,n,r,var,x) [

Function() SumFunc(k,from,to,summand, sum, predicate );
Function() SumFunc(k,from,to,summand, sum);
HoldArg(SumFunc,predicate);
HoldArg(SumFunc,sum);
HoldArg(SumFunc,summand);

// Difference code does not work
SumFunc(_sumvar,sumfrom_IsInteger,_sumto,_sumbody,_sum) <--
[
	// Take the given answer and create 2 rules, one for an exact match
	// for sumfrom, and one which will catch sums starting at a different
	// index and subtract off the difference

	`(40 # Sum(@sumvar,@sumfrom,@sumto,@sumbody )	<-- Eval(@sum) );
	`(41 # Sum(@sumvar,p_IsInteger,@sumto,@sumbody)_(p > @sumfrom) 
	     <-- 
	     [
	          Local(sub);
		  (sub := Eval(UnList({Sum,sumvar'arg,@sumfrom,p-1,sumbody'arg})));
		  Simplify(Eval(@sum) - sub );
             ]);
];

SumFunc(_sumvar,sumfrom_IsInteger,_sumto,_sumbody,_sum,_condition) <--
[

	`(40 # Sum(@sumvar,@sumfrom,@sumto,@sumbody)_(@condition)    <-- Eval(@sum) );
	`(41 # Sum(@sumvar,p_IsInteger,@sumto,@sumbody )_(@condition And p > @sumfrom) 
	     <-- 
	     [
	          Local(sub);
		  `(sub := Eval(UnList({Sum,sumvar'arg,@sumfrom,p-1,sumbody'arg})));
		  Simplify(Eval(@sum) - sub );
             ]);
];

// Some type of canonical form is needed so that these match when
// given in a different order, like x^k/k! vs. (1/k!)*x^k 
// works !
SumFunc(_k,1,_n,_c + _d,
  Eval(UnList({Sum,sumvar'arg,1,n,c})) + 
  Eval(UnList({Sum,sumvar'arg,1,n,d})) 
);
SumFunc(_k,1,_n,_c*_expr,Eval(c*UnList({Sum,sumvar'arg,1,n,expr})), IsFreeOf(k,c) );
SumFunc(_k,1,_n,_expr/_c,Eval(UnList({Sum,sumvar'arg,1,n,expr})/c), IsFreeOf(k,c) );

// this only works when the index=1
// If the limit of the general term is not zero, then the series diverges
// We need something like IsUndefined(term), because this croaks when limit return Undefined
//SumFunc(_k,1,Infinity,_expr,Infinity,Eval(Abs(UnList({Limit,sumvar'arg,Infinity,expr})) > 0));
SumFunc(_k,1,Infinity,1/_k,Infinity);

SumFunc(_k,1,_n,_c,c*n,IsFreeOf(k,c) );
SumFunc(_k,1,_n,_k, n*(n+1)/2 );
//SumFunc(_k,1,_n,_k^2, n*(n+1)*(2*n+1)/6 );
//SumFunc(_k,1,_n,_k^3, (n*(n+1))^2 / 4 );
SumFunc(_k,1,_n,_k^_p,(Bernoulli(p+1,n+1) - Bernoulli(p+1))/(p+1), IsInteger(p) );
SumFunc(_k,1,_n,2*_k-1, n^2 );
SumFunc(_k,1,_n,HarmonicNumber(_k),(n+1)*HarmonicNumber(n) - n );

// Geometric series! The simplest of them all ;-)
SumFunc(_k,0,_n,(r_IsFreeOf(k))^(_k), (1-r^(n+1))/(1-r) );

// Infinite Series
// this allows Zeta a complex argument, which is not supported yet
SumFunc(_k,1,Infinity,1/(_k^_d), Zeta(d), IsFreeOf(k,d) );
SumFunc(_k,1,Infinity,_k^(-_d), Zeta(d), IsFreeOf(k,d) );

SumFunc(_k,0,Infinity,_x^(2*_k+1)/(2*_k+1)!,Sinh(x) );
SumFunc(_k,0,Infinity,(-1)^k*_x^(2*_k+1)/(2*_k+1)!,Sin(x) );
SumFunc(_k,0,Infinity,_x^(2*_k)/(2*_k)!,Cosh(x) );
SumFunc(_k,0,Infinity,(-1)^k*_x^(2*_k)/(2*_k)!,Cos(x) );
SumFunc(_k,0,Infinity,_x^(2*_k+1)/(2*_k+1),ArcTanh(x) );
SumFunc(_k,0,Infinity,1/(_k)!,Exp(1) );
SumFunc(_k,0,Infinity,_x^_k/(_k)!,Exp(x) );
40 # Sum(_var,_from,Infinity,_expr)_( `(Limit(@var,Infinity)(@expr)) = Infinity) <-- Infinity;

SumFunc(_k,1,Infinity,1/Bin(2*_k,_k), (2*Pi*Sqrt(3)+9)/27 );
SumFunc(_k,1,Infinity,1/(_k*Bin(2*_k,_k)), (Pi*Sqrt(3))/9 );
SumFunc(_k,1,Infinity,1/(_k^2*Bin(2*_k,_k)), Zeta(2)/3 );
SumFunc(_k,1,Infinity,1/(_k^3*Bin(2*_k,_k)), 17*Zeta(4)/36 );
SumFunc(_k,1,Infinity,(-1)^(_k-1)/_k, Ln(2) );

];