Menu

#5187 rootscontract(%i*sqrt(4*%i-2)) is wrong

None
open
nobody
5
3 days ago
3 days ago
No
(%i1) rootscontract(%i*sqrt(4*%i-2));
(%o1) sqrt(2-4*%i)

The correct answer is -sqrt(2-4*%i).

Discussion

  • David Scherfgen

    David Scherfgen - 3 days ago

    This causes a wrong result in integrate(asin(2+%i*x),x,-1,1).

     
  • Stavros Macrakis

    Are you sure?

    %i*sqrt(4*%i-2)) 
       == sqrt(-1)*sqrt(4*%i-2))
       contract =>
       sqrt( (-1)*(4*%i-2) )
       ==
       sqrt(2-4*%i)
    
     
    • David Scherfgen

      David Scherfgen - 3 days ago

      By that logic:

      %i * %i = sqrt(-1) * sqrt(-1)
      contract => sqrt((-1) * (-1)) = sqrt(1) = 1
      

      I have a patch: about ten lines in rtcon / rtc-fixitup (src/comm2.lisp). It leaves contraction of an actual product of roots alone, and drops two rewrites that aren't contractions at all:

      • %i, which rt-separ collects as the root (-1)^(1/2). As a radicand, -1 has argument %pi — the whole (-%pi, %pi] budget — so any companion with a positive argument wraps. And nothing is gained: for a positive companion the result is just %i*sqrt(2) respelled as sqrt(-2).
      • a lone root, rewritten R^(N/D) -> (R^N)^(1/D). No product to contract, and it wraps once N*arg(R) leaves (-%pi, %pi]:
      (%i1) rootscontract(2*(-1)^(3/4));
      (%o1)                       2*(-1)^(1/4)
      

      Besides the limit bug it also fixes algsys([x^4+1],[x]), which gave sqrt(-%i) as its fourth root — that is -(-1)^(1/4)*%i over again, three distinct roots for a quartic — to_poly_solve(x^8=1,x) losing %i/sqrt(2)-1/sqrt(2), and the wrong branch asserted by problems 308 and 309 of rtest_limit_extra.mac. Rest of the suite unchanged.

      It does not make rootscontract sound, though: sqrt(-2+4*%i)*sqrt(-3+%i) still comes back negated.

       
  • Stavros Macrakis

    Yes, rootscontract in general doesn't preserve the choice of sign. I think that's inherent in its functionality, isn't it?

     
    • David Scherfgen

      David Scherfgen - 3 days ago

      For a genuine product of roots, yes - I'd accept that as inherent, and the patch leaves it exactly as it is.

      But neither of the two cases it removes is a contraction:

      • %i is not a root in the expression - rt-separ "invents" (-1)^(1/2) and folds it in.
      • a lone root, R^(N/D) -> (R^N)^(1/D), has no product to contract.

      Nothing gets combined in either case, so there is no functionality to preserve — only the sign risk, taken for free.

      And it isn't only about the spelling of a user-visible answer, since rootscontract is called internally: algsys([x^4+1],[x]) returns sqrt(-%i) as its fourth root, which is -(-1)^(1/4)*%i listed again, for a second time - a quartic with three distinct roots - and to_poly_solve(x^8=1,x) drops %i/sqrt(2)-1/sqrt(2). So either rootscontract gets stricter where it gains nothing, or its internal callers shouldn't be using it.

       
  • Stavros Macrakis

    Consider

    rootscontract(sqrt(a)*sqrt(4*%i-2))
    

    Should that produce a different result when you substitute -1 for a?

     
    • David Scherfgen

      David Scherfgen - 3 days ago

      True, they differ under the proposed patch:

      subst(a=-1, rootscontract(sqrt(a)*sqrt(4*%i-2)))  =>  sqrt(2-4*%i)     =  1.7989 - 1.1118*%i
      rootscontract(%i*sqrt(4*%i-2))                    =>  %i*sqrt(4*%i-2)  = -1.7989 + 1.1118*%i
      

      Without the patch, both give sqrt(2-4*%i). But that agreement is agreement on a wrong value - the expression being rewritten, %i*sqrt(4*%i-2), is -sqrt(2-4*%i). So the old behaviour preserved the value along neither route; it failed identically along both. The patch leaves the first route as it was and makes the second one correct.

      Getting the two routes to agree and be right needs a stricter rootscontract, which I'd suggest adding as an opt-in rather than changing the default. Concretely, a new option variable

      rootscontract_safe (boolean)
      default: false
      

      kept separate from rootsconmode, since that one governs which denominators combine - an orthogonal axis, and you'd want safe crossed with all three of its values. At false everything behaves exactly as it does now. At true, rootscontract performs only contractions it can prove preserve the value, and leaves the rest uncontracted; it never contracts more than the default does.

      The reason for defaulting to false is that the internal callers in limit.lisp and algsys.lisp can bind it to true, so limit and algsys stop handing out wrong values, while a direct rootscontract(...) by the user behaves as it always has. Compatibility is a concern, I guess.

      On your previous example: 4*%i-2 is not real, so with rootscontract_safe:true nothing is provable and both routes are refused - both come out %i*sqrt(4*%i-2). Under assume(a>0) both contract, and a=-1 no longer arises.

       

Log in to post a comment.