f2c04ccdad 2021-03-03 trnsz@pobox.c: %%%%%%%%%%%%%%%%%%%%%
f2c04ccdad 2021-03-03 trnsz@pobox.c: % PROGRAMMING
f2c04ccdad 2021-03-03 trnsz@pobox.c: %%%%%%%%%%%%%%%%%%%%%
f2c04ccdad 2021-03-03 trnsz@pobox.c:
f2c04ccdad 2021-03-03 trnsz@pobox.c: % Define number to factorize
f2c04ccdad 2021-03-03 trnsz@pobox.c: x:=42;
f2c04ccdad 2021-03-03 trnsz@pobox.c:
f2c04ccdad 2021-03-03 trnsz@pobox.c: % Factorize x and write out each individual
f2c04ccdad 2021-03-03 trnsz@pobox.c: % factor
f2c04ccdad 2021-03-03 trnsz@pobox.c: factors:=factorize(fix(x))$
f2c04ccdad 2021-03-03 trnsz@pobox.c: x:=0$
f2c04ccdad 2021-03-03 trnsz@pobox.c: for i:=1:length(factors) do begin
f2c04ccdad 2021-03-03 trnsz@pobox.c: q:=part(factors,i);
f2c04ccdad 2021-03-03 trnsz@pobox.c: for j:=1:part(q,2) do begin
f2c04ccdad 2021-03-03 trnsz@pobox.c: x:=x+1;
f2c04ccdad 2021-03-03 trnsz@pobox.c: write "factor ", x, ": ", part(q,1);
f2c04ccdad 2021-03-03 trnsz@pobox.c: end;
f2c04ccdad 2021-03-03 trnsz@pobox.c: end;
f2c04ccdad 2021-03-03 trnsz@pobox.c:
f2c04ccdad 2021-03-03 trnsz@pobox.c: % Procedure to calculate Legendre polynomial
f2c04ccdad 2021-03-03 trnsz@pobox.c: % using recursion
f2c04ccdad 2021-03-03 trnsz@pobox.c: procedure p(n,x);
f2c04ccdad 2021-03-03 trnsz@pobox.c: if n<0 then rederr "Invalid argument to p(n,x)"
f2c04ccdad 2021-03-03 trnsz@pobox.c: else if n=0 then 1
f2c04ccdad 2021-03-03 trnsz@pobox.c: else if n=1 then x
f2c04ccdad 2021-03-03 trnsz@pobox.c: else ((2*n-1)*x*p(n-1,x)-(n-1)*p(n-2,x))/n$
f2c04ccdad 2021-03-03 trnsz@pobox.c:
f2c04ccdad 2021-03-03 trnsz@pobox.c: % Enable fancy output
f2c04ccdad 2021-03-03 trnsz@pobox.c: fancy_output$
f2c04ccdad 2021-03-03 trnsz@pobox.c:
f2c04ccdad 2021-03-03 trnsz@pobox.c: % Calculate p(2,w)
f2c04ccdad 2021-03-03 trnsz@pobox.c: write "P(2,w) = ", p(2,w);
f2c04ccdad 2021-03-03 trnsz@pobox.c:
f2c04ccdad 2021-03-03 trnsz@pobox.c: % Incidentally, p(n,x) can be calculated more
f2c04ccdad 2021-03-03 trnsz@pobox.c: % efficiently as follows
f2c04ccdad 2021-03-03 trnsz@pobox.c: procedure p(n,x);
f2c04ccdad 2021-03-03 trnsz@pobox.c: sub(y=0,df(1/(y^2-2*x*y+1)^(1/2),y,n))/(for i:=1:n product i)$
f2c04ccdad 2021-03-03 trnsz@pobox.c:
f2c04ccdad 2021-03-03 trnsz@pobox.c: write "P(3,w) = ", p(3,w);