critical_points(expr,var):= block([numer:true,f:diff(expr,var)], realroots(num(f)*denom(f),10^-7)); function_max(expr,var,a,b,critical_points):= block([y:max(at(expr,var=a),at(expr,var=b)) ], for v in critical_points do (v:at(x,v), if a y then y:at(expr,var=v)), y); upper_sum(expr,var,partition):= block([crit:critical_points(expr,var)], partition:sort(partition), sum(function_max(expr,var,partition[i],partition[i+1],crit)* (partition[i+1]-partition[i]),i,1,length(partition)-1)); /* this will be part of maxima in a little while .. */ graph_riemann_sum(up,expr,var,partition):= block([lis,numer:true,crit:critical_points(expr,var),n:length(partition),xx, interval,display2d:false,sgn,min:100000,max:-10000000], sgn: if up= upper then 1 else -1, partition:sort(partition), interval:partition[n]-partition[1], lis:[], for i:0 thru 50 do (xx:partition[1]+i*interval/50, lis:cons(val:at(expr,var=xx),lis), lis:cons(xx,lis), if val < min then min:val, if val > max then max:val), sprint("{plot2d ", "{xrange " ,partition[1]-.5,partition[n]+.5,"} {yrange" , min-.5,max+.5, "} {xversusy "), tcl_output(lis,1), tcl_output(lis,2), sprint( "}"), sprint(" {xversusy {"), for i:1 thru n-1 do( sprint(partition[i],partition[i],partition[i+1],partition[i+1], partition[i])), sprint(" } { "), for i:1 thru n-1 do(xx:sgn*function_max(sgn*expr,var,partition[i],partition[i+1],crit), sprint(0.0,xx,xx,0.0,0.0)), sprint("}}}"), "" ); upper_and_lower_sums(expr,var,partition):= block([up:upper_sum(expr,var,partition), low:-upper_sum(-expr,var,partition)], [up,low,up-low]); make_partition(a,b,n):=makelist(a+(b-a)*i/n,i,0,n);