[Mpc-discuss] mpc_pow(big real, big imag) gives incorrect output at low precision

Simon Tatham Simon.Tatham at arm.com
Tue Oct 14 12:27:04 CEST 2014


The attached source file uses mpc_pow() to raise 2^1000 to the power (2^1000 i), in a range of working precisions from 100 bits to 2000 bits.

As I understand MPC, the output value should be computed as if the true exact answer was first found, and then rounded to the nearest representable value in the output precision. Therefore, if I print the first few significant digits of the result, I should expect them to be consistent across all the working precisions in which I attempted the calculation.

In fact, when I run the program, I see the following output:

$ gcc -o complex-pow-bug complex-pow-bug.c -lmpc -lmpfr -lgmp && ./complex-pow-bug
prec =  100 : -0.65998018613219828377 + -0.75128300520703175837 i
prec =  200 : -0.37810342220239686510 + -0.92576336183650950389 i
prec =  300 : -0.88265426875960838977 + 0.47002280991505188965 i
prec =  400 : 0.75398492012629596058 + 0.65689172640713256324 i
prec =  500 : -0.23471860577371806070 + 0.97206336012815642318 i
prec =  600 : 0.03552424414330536381 + 0.99936881484167137011 i
prec =  700 : 0.00525482109750118460 + 0.99998619333230459180 i
prec =  800 : 0.88177125100583345985 + 0.47167728469750103562 i
prec =  900 : 0.43001394983546122677 + 0.90282224327212133414 i
prec = 1000 : 0.99530388965077993064 + 0.09679962420396107553 i
prec = 1100 : 0.99530388965077993064 + 0.09679962420396107553 i
prec = 1200 : 0.99530388965077993064 + 0.09679962420396107553 i
prec = 1300 : 0.99530388965077993064 + 0.09679962420396107553 i
prec = 1400 : 0.99530388965077993064 + 0.09679962420396107553 i
prec = 1500 : 0.99530388965077993064 + 0.09679962420396107553 i
prec = 1600 : 0.99530388965077993064 + 0.09679962420396107553 i
prec = 1700 : 0.99530388965077993064 + 0.09679962420396107553 i
prec = 1800 : 0.99530388965077993064 + 0.09679962420396107553 i
prec = 1900 : 0.99530388965077993064 + 0.09679962420396107553 i
prec = 2000 : 0.99530388965077993064 + 0.09679962420396107553 i

in which the output varies wildly until the precision reaches 1000 bits, and then settles down to what I presume is the actually correct answer.

This particular case of complex pow() requires the computation of (2^1000 log 2^1000) to very high precision, so that it can be range-reduced for sin and cos and still have enough precision left for the output. From the fact that the output stabilised at 1000 bits of precision, I'm guessing that the cause of this behaviour is that the logarithm is computed in only the input precision, rather than using higher precision internally.

I observed this on an Ubuntu 12.04 x86_64 machine, with the following compiler and libraries:
 - gcc from the Ubuntu 'gcc' package, version 4:4.6.3-1ubuntu5. 'gcc -v' reports
     gcc version 4.6.3 (Ubuntu/Linaro 4.6.3-1ubuntu5)
 - MPFR from the Ubuntu 'libmpfr4' package, version 3.1.0-3ubuntu2
 - GMP from the Ubuntu 'libgmp10' package, version 2:5.0.2+dfsg-2ubuntu1
 - 'uname -a' reports
     Linux e103988-lin 3.2.0-70-generic #105-Ubuntu SMP Wed Sep 24 19:49:16 UTC 2014 x86_64 x86_64 x86_64 GNU/Linux

and I saw identical output from both the Ubuntu-packaged version of MPC (package 'libmpc2' version 0.9-4) and from the up-to-date MPC 1.0.2 which I compiled by hand from the source tarball.

Cheers,
Simon

-- IMPORTANT NOTICE: The contents of this email and any attachments are confidential and may also be privileged. If you are not the intended recipient, please notify the sender immediately and do not disclose the contents to any other person, use it for any purpose, or store or copy the information in any medium.  Thank you.

ARM Limited, Registered office 110 Fulbourn Road, Cambridge CB1 9NJ, Registered in England & Wales, Company No:  2557590
ARM Holdings plc, Registered office 110 Fulbourn Road, Cambridge CB1 9NJ, Registered in England & Wales, Company No:  2548782
-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: complex-pow-bug.c
URL: <http://lists.gforge.inria.fr/pipermail/mpc-discuss/attachments/20141014/a08e9ddc/attachment.c>


More information about the Mpc-discuss mailing list