simon on Wed, 18 Dec 2013 11:47:52 +0100


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Re: Polresultant bug ?



Hi,

This is not a problem of polresultant, but of precision in your computations.

Indeed, the coefficients of your annulator_for_c have about 80 digits.
In order to have a coherent result, you have to do your computations with a number of digits significatively larger than 80.
Try for example 2*80 = 160.

Denis SIMON.



On Wed, 18 Dec 2013, Ewan Delanoy wrote:



  Hello all,

  it seems I found an example where according to GP, the resultant of two
expressions that evaluate
to zero gives an expression that evaluate to a nonzero (and in fact, huge)
value.
  Below I paste both the source code and the output I got.
  Thanks in advance for any help,

  Ewan


/****************************************************************************
********** */
/*The source code */
/****************************************************************************
********** */

annulator_for_a=256*a^12 - 363776*a^10 + 193847136*a^8 - 46812817968*a^6 +
4719167653697*a^4 - 114005664356542*a^2 + 1412780238606961
floating_a=polroots(annulator_for_a)[1]
annulator_for_b=b^6 - 108*b^4 + 2684*b^2 - 59536
floating_b=real(polroots(annulator_for_b)[2])
floating_c=floating_a-floating_b
term1=subst(annulator_for_a,a,b+c)
should_be_zero1=subst(subst(term1,b,floating_b),c,floating_c)
annulator_for_c=polresultant(term1,annulator_for_b,b)
should_be_zero2=subst(annulator_for_c,c,floating_c)

/****************************************************************************
********** */
/*The output  */
/****************************************************************************
********** */

MacBook-Pro-de-Ewan:~ ewandelanoy$ gp
                  GP/PARI CALCULATOR Version 2.5.3 (released)
          i386 running darwin (x86-64/GMP-5.0.4 kernel) 64-bit version
compiled: Jan 24 2013, gcc-4.2.1 (Based on Apple Inc. build 5658) (LLVM
build 2336.11.00)
                 (readline v6.2 enabled, extended help enabled)

                     Copyright (C) 2000-2011 The PARI Group

PARI/GP is free software, covered by the GNU General Public License, and
comes
WITHOUT ANY WARRANTY WHATSOEVER.

Type ? for help, \q to quit.
Type ?12 for how to get moral (and possibly technical) support.

parisize = 8000000, primelimit = 500509
? annulator_for_a=256*a^12 - 363776*a^10 + 193847136*a^8 - 46812817968*a^6 +
4719167653697*a^4 - 114005664356542*a^2 + 1412780238606961
%1 = 256*a^12 - 363776*a^10 + 193847136*a^8 - 46812817968*a^6 +
4719167653697*a^4 - 114005664356542*a^2 + 1412780238606961
? floating_a=polroots(annulator_for_a)[1]
%2 = 16.507639886376369964102375152715808868 -
1.1732991791154715940767062706818971012*I
? annulator_for_b=b^6 - 108*b^4 + 2684*b^2 - 59536
%3 = b^6 - 108*b^4 + 2684*b^2 - 59536
? floating_b=real(polroots(annulator_for_b)[2])
%4 = 9.1973372100322392595456421736424968201
? floating_c=floating_a-floating_b
%5 = 7.3103026763441307045567329790733120478 -
1.1732991791154715940767062706818971012*I
? term1=subst(annulator_for_a,a,b+c)
%6 = 256*b^12 + 3072*c*b^11 + (16896*c^2 - 363776)*b^10 + (56320*c^3 -
3637760*c)*b^9 + (126720*c^4 - 16369920*c^2 + 193847136)*b^8 + (202752*c^5 -
43653120*c^3 + 1550777088*c)*b^7 + (236544*c^6 - 76392960*c^4 +
5427719808*c^2 - 46812817968)*b^6 + (202752*c^7 - 91671552*c^5 +
10855439616*c^3 - 280876907808*c)*b^5 + (126720*c^8 - 76392960*c^6 +
13569299520*c^4 - 702192269520*c^2 + 4719167653697)*b^4 + (56320*c^9 -
43653120*c^7 + 10855439616*c^5 - 936256359360*c^3 + 18876670614788*c)*b^3 +
(16896*c^10 - 16369920*c^8 + 5427719808*c^6 - 702192269520*c^4 +
28315005922182*c^2 - 114005664356542)*b^2 + (3072*c^11 - 3637760*c^9 +
1550777088*c^7 - 280876907808*c^5 + 18876670614788*c^3 -
228011328713084*c)*b + (256*c^12 - 363776*c^10 + 193847136*c^8 -
46812817968*c^6 + 4719167653697*c^4 - 114005664356542*c^2 +
1412780238606961)
? should_be_zero1=subst(subst(term1,b,floating_b),c,floating_c)
%7 = 9.556391613824554411 E-22 - 3.105697067803841053 E-22*I
? annulator_for_c=polresultant(term1,annulator_for_b,b)
%8 = 281474976710656*c^72 - 2764647221252063232*c^70 +
12881090734407523762176*c^68 - 37900879189576872787705856*c^66 +
79099815954946491364482744320*c^64 - 124684759166630772970798278770688*c^62
+ 154348991980763936851992593913872384*c^60 -
154045742043145589406276188287109431296*c^58 +
126278018863764650417996510658035946881024*c^56 -
86188295361347830927179212093084919761207296*c^54 +
49478930044049671150743174794661019965550755840*c^52 -
24074514186962292770585768259851360789778778816512*c^50 +
9984746464975367599334473741516071684995601427070976*c^48 -
3544601773895860008395886389523866780061142280317173760*c^46 +
1080184220852605057051843703838564046424092897453830832128*c^44 -
283063744081987480169735131103108552441984386863102852857856*c^42 +
63829717334668411467356121108379367011430103612372083744374784*c^40 -
12380931236742739255452291376833192852560602832000266561304526848*c^38 +
2062669465587519366960798626747033680113649744703395573266979635200*c^36 -
294342896040799822278685004581163889813837967268955216247392575873024*c^34 +35822170668621222716189038834945741067330091385125921979038859956748800*c^3
2 -3694599094972804492903490053128068674286670925141467454235954938531409408*c
^30 +319960240851023753657282679654500912353276531225489928628714114061458570560
*c^28 -229523124348646894674693749507943423678604971743540760886488234477385582240
32*c^26 +133540122622880453154192069395728174480573210717450941020933452226076490321
5361*c^24 -608091485575824049928810016694298694113579523852570953890936245732723388993
76068*c^22 +201756641433254568341552129401331699270034937209354126404620590151755348180
8626562*c^20 -396135443048892631637411626519338142317965576647394035113621469768341039106
76919316*c^18 -830734957146435644175388145522917471286550249793111950753257626172855716538
76743441*c^16 +332312675119643619106768952565308776883182037008020918543610524912665691806
88206374776*c^14 -990437872605966242896879211543870708900639163340555827418555490651344092124
153842082212*c^12 +932912335531598587420267492196329214620107896529499931168264848122095697606
4727122600216*c^10 +175291809114698568461616198975915940771681390157547657579559269812558585841
725448099842799*c^8 -551646020293422160604256097389154281421501271527752864787625977318114545794
3877549960013748*c^6 +418288683874030712474395958976663314737499086339637924696074115753271457709
62254753636273922*c^4 +839876356297955607038827371458654374227011532016895487424475243901159477944
48018764725358492*c^2 +503773428948680856170140977878135888394567173775246602443797048650457162370
78210026453561601
? should_be_zero2=subst(annulator_for_c,c,floating_c)
%9 = 2.2768496195346216922 E60 - 5.226873626511962736 E59*I