[GiNaC-devel] Simplifying powers

[Alexei Sheplyakov]

Hi again,

> +	ex e1 = pow(5*pow(3*a*b*x*y*p*q,2),7*half*c).expand();
> +	ex e2 = pow(p,7*c)*pow(q,7*c)*pow(pow(a*b*x*y,2),numeric(7,2)*c)*pow(45,numeric(7,2)*c);

The idea looks good...

> diff --git a/ginac/power.cpp b/ginac/power.cpp
> index 9a8f7d5..d01a88f 100644
> --- a/ginac/power.cpp
> +++ b/ginac/power.cpp
> @@ -824,11 +824,38 @@ ex power::expand(unsigned options) const
>  	
>  	if (!is_exactly_a<numeric>(expanded_exponent) ||
>  		!ex_to<numeric>(expanded_exponent).is_integer()) {
> -		if (are_ex_trivially_equal(basis,expanded_basis) && are_ex_trivially_equal(exponent,expanded_exponent)) {
> -			return this->hold();
> -		} else {
> +		if (is_exactly_a<mul>(expanded_basis) && expanded_exponent.info(info_flags::real)) {

... however the implementation is a bit inefficient. Even worse -- it makes
power::expand() slower even if the above transformation is not applicable.
Consider the following expression:

((1+x)^(2*eps)*(1+2*x)^(1+2*eps)*...(1+100*x)^(1+100*eps))^(4-2*eps)

> +			const mul &m = ex_to<mul>(expanded_basis);
> +			exvector prodseq, powseq;
> +			prodseq.reserve(m.seq.size() + 1);
> +			powseq.reserve(m.seq.size() + 1);
> +			epvector::const_iterator last = m.seq.end();
> +			epvector::const_iterator cit = m.seq.begin();
> +			bool possign=true;
> +
> +			while (cit!=last) {
> +				ex e=m.recombine_pair_to_ex(*cit);
> +				if (e.info(info_flags::positive))
> +					prodseq.push_back(power(e, expanded_exponent));
> +				else if (e.info(info_flags::negative)) {
> +					prodseq.push_back(power(-e, expanded_exponent));
> +					possign = !possign;

So we check every term here (which might be quite expensive)...

> +				} else
> +					powseq.push_back(m.recombine_pair_to_ex(*cit));

This conversion is useless (mul(powseq) converts the exvector back to epvector).
Perhaps powseq should be a epvector.

> +				++cit;
> +			}
> +
> +			if (m.overall_coeff.info(info_flags::positive))
> +				prodseq.push_back(power(m.overall_coeff, expanded_exponent));
> +			else if (m.overall_coeff.info(info_flags::negative)) {
> +				prodseq.push_back(power(-m.overall_coeff, expanded_exponent));
> +				possign = !possign;
> +			} else
> +				powseq.push_back(m.overall_coeff);
> +
> +			return (new mul(prodseq))->setflag(status_flags::dynallocated)*(new power((possign? _ex1 : _ex_1)*mul(powseq), expanded_exponent))->setflag(status_flags::dynallocated);

Eventually the loop constructs the same expression, but it needs to be
evaluated once again.

> +		} else
>  			return (new power(expanded_basis,expanded_exponent))->setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0));
> -		}
>  	}

Therefore I think the patch should be improved. In particular,

- info_flags::{positive,negative} should be cached (add yet another
  status_flags and maintain them)
- the trivial case should be handled more carefully (i.e. the code should
  detect that the transformation is not applicable as early as possible,
  and return the original expression).

Best regards,
	Alexei