Annotation of OpenXM_contrib/pari-2.2/src/basemath/subcyclo.c, Revision 1.1
1.1 ! noro 1: /* $Id: subcyclo.c,v 1.21 2002/07/04 12:25:04 bill Exp $
! 2:
! 3: Copyright (C) 2000 The PARI group.
! 4:
! 5: This file is part of the PARI/GP package.
! 6:
! 7: PARI/GP is free software; you can redistribute it and/or modify it under the
! 8: terms of the GNU General Public License as published by the Free Software
! 9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
! 10: ANY WARRANTY WHATSOEVER.
! 11:
! 12: Check the License for details. You should have received a copy of it, along
! 13: with the package; see the file 'COPYING'. If not, write to the Free Software
! 14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
! 15:
! 16: #include "pari.h"
! 17:
! 18: extern GEN vandermondeinversemod(GEN L, GEN T, GEN den, GEN mod);
! 19: extern GEN supnorm(GEN L, long prec);
! 20:
! 21: /*************************************************************************/
! 22: /** **/
! 23: /** Routines for handling subgroups of (Z/nZ)^* **/
! 24: /** without requiring discrete logarithms. **/
! 25: /** **/
! 26: /*************************************************************************/
! 27:
! 28: /* Subgroups are [gen,ord,bits] where
! 29: * gen is a vecsmall of generators
! 30: * ord is theirs relative orders
! 31: * bits is a bit vector of the elements, of length(n).
! 32: */
! 33:
! 34: /*The algorithm is similar to testpermutation*/
! 35: void
! 36: znstar_partial_coset_func(long n, GEN H, void (*func)(void *data,long c)
! 37: , void *data, long d, long c)
! 38: {
! 39: GEN gen = (GEN) H[1];
! 40: GEN ord = (GEN) H[2];
! 41: GEN cache = vecsmall_const(d,c);
! 42: long i, j, card = 1;
! 43:
! 44: (*func)(data,c);
! 45: for (i = 1; i <= d; i++) card *= ord[i];
! 46: for(i=1; i<card; i++)
! 47: {
! 48: long k, m = i;
! 49: for(j=1; j<d && m%ord[j]==0 ;j++) m /= ord[j];
! 50: cache[j] = mulssmod(cache[j],gen[j],n);
! 51: for (k=1; k<j; k++) cache[k] = cache[j];
! 52: (*func)(data, cache[j]);
! 53: }
! 54: }
! 55:
! 56: void
! 57: znstar_coset_func(long n, GEN H, void (*func)(void *data,long c)
! 58: , void *data, long c)
! 59: {
! 60: znstar_partial_coset_func(n, H, func,data, lg(H[1])-1, c);
! 61: }
! 62:
! 63: /*Add the element of the bitvec of the coset c modulo the subgroup of H
! 64: * generated by the first d generators to the bitvec bits.*/
! 65:
! 66: void
! 67: znstar_partial_coset_bits_inplace(long n, GEN H, GEN bits, long d, long c)
! 68: {
! 69: gpmem_t ltop=avma;
! 70: znstar_partial_coset_func(n,H, (void (*)(void *,long)) &bitvec_set,
! 71: (void *) bits, d, c);
! 72: avma=ltop;
! 73: }
! 74:
! 75: void
! 76: znstar_coset_bits_inplace(long n, GEN H, GEN bits, long c)
! 77: {
! 78: znstar_partial_coset_bits_inplace(n, H, bits, lg(H[1])-1, c);
! 79: }
! 80:
! 81: GEN
! 82: znstar_partial_coset_bits(long n, GEN H, long d, long c)
! 83: {
! 84: GEN bits=bitvec_alloc(n);
! 85: znstar_partial_coset_bits_inplace(n,H,bits,d,c);
! 86: return bits;
! 87: }
! 88:
! 89: /*Compute the bitvec of the elements of the subgroup of H generated by the
! 90: * first d generators.
! 91: */
! 92:
! 93: GEN
! 94: znstar_coset_bits(long n, GEN H, long c)
! 95: {
! 96: return znstar_partial_coset_bits(n, H, lg(H[1])-1, c);
! 97: }
! 98:
! 99: /*Compute the bitvec of the elements of the subgroup of H generated by the
! 100: * first d generators.*/
! 101:
! 102: GEN
! 103: znstar_partial_bits(long n, GEN H, long d)
! 104: {
! 105: return znstar_partial_coset_bits(n, H, d, 1);
! 106: }
! 107: /*Compute the bitvec of the elements of H.*/
! 108:
! 109: GEN
! 110: znstar_bits(long n, GEN H)
! 111: {
! 112: return znstar_partial_bits(n,H,lg(H[1])-1);
! 113: }
! 114:
! 115: /*Compute the subgroup of (Z/nZ)^* generated by the elements of
! 116: * the vecsmall V.
! 117: */
! 118:
! 119: GEN
! 120: znstar_generate(long n, GEN V)
! 121: {
! 122: gpmem_t ltop=avma;
! 123: GEN res=cgetg(4,t_VEC);
! 124: GEN gen=cgetg(lg(V),t_VECSMALL);
! 125: GEN ord=cgetg(lg(V),t_VECSMALL);
! 126: GEN bits;
! 127: long i,r=0;
! 128: res[1]=(long)gen;
! 129: res[2]=(long)ord;
! 130: bits=znstar_partial_bits(n,res,r);
! 131: for(i=1;i<lg(V);i++)
! 132: {
! 133: long v=V[i];
! 134: long g=v;
! 135: long o=0;
! 136: while(!bitvec_test(bits,g))
! 137: {
! 138: g=mulssmod(g,v,n);
! 139: o++;
! 140: }
! 141: if (o)
! 142: {
! 143: r++;
! 144: gen[r]=v;
! 145: ord[r]=o+1;
! 146: cgiv(bits);
! 147: bits=znstar_partial_bits(n,res,r);
! 148: }
! 149: }
! 150: setlg(gen,r+1);
! 151: setlg(ord,r+1);
! 152: res[3]=(long)bits;
! 153: return gerepilecopy(ltop,res);
! 154: }
! 155:
! 156: /* Return the lists of element of H.
! 157: * This can be implemented with znstar_coset_func instead.
! 158: */
! 159:
! 160: GEN
! 161: znstar_elts(long n, GEN H)
! 162: {
! 163: long card=group_order(H);
! 164: GEN gen=(GEN)H[1], ord=(GEN)H[2];
! 165: GEN sg = cgetg(1 + card, t_VECSMALL);
! 166: long k, j, l;
! 167: sg[1] = 1;
! 168: for (j = 1, l = 1; j < lg(gen); j++)
! 169: {
! 170: int c = l * (ord[j] - 1);
! 171: for (k = 1; k <= c; k++) /* I like it */
! 172: sg[++l] = mulssmod(sg[k], gen[j], n);
! 173: }
! 174: vecsmall_sort(sg);
! 175: return sg;
! 176: }
! 177:
! 178: /* Take a znstar H and n dividing the modulus of H.
! 179: * Output H reduced to modulus n */
! 180:
! 181: GEN
! 182: znstar_reduce_modulus(GEN H, long n)
! 183: {
! 184: gpmem_t ltop=avma;
! 185: GEN gen=cgetg(lg(H[1]),t_VECSMALL);
! 186: long i;
! 187: for(i=1; i < lg(gen); i++)
! 188: gen[i] = mael(H,1,i)%n;
! 189: return gerepileupto(ltop, znstar_generate(n,gen));
! 190: }
! 191: /*Compute conductor of H*/
! 192:
! 193: long znstar_conductor(long n, GEN H)
! 194: {
! 195: gpmem_t ltop=avma;
! 196: int i,j;
! 197: GEN F;
! 198: long cnd=n;
! 199: F = decomp_small(n);
! 200: for(i=lg((GEN)F[1])-1;i>0;i--)
! 201: {
! 202: long p=coeff(F,i,1);
! 203: long e=coeff(F,i,2);
! 204: long q=n;
! 205: if (DEBUGLEVEL>=4)
! 206: fprintferr("SubCyclo: testing %ld^%ld\n",p,e);
! 207: for ( ; e>=1; e--)
! 208: {
! 209: long z = 1;
! 210: q /= p;
! 211: for (j = 1; j < p; j++)
! 212: {
! 213: z += q;
! 214: if (!bitvec_test((GEN) H[3],z) && cgcd(z,n)==1)
! 215: break;
! 216: }
! 217: if ( j < p )
! 218: {
! 219: if (DEBUGLEVEL>=4)
! 220: fprintferr("SubCyclo: %ld not found\n",z);
! 221: break;
! 222: }
! 223: cnd /= p;
! 224: if (DEBUGLEVEL>=4)
! 225: fprintferr("SubCyclo: new conductor:%ld\n",cnd);
! 226: }
! 227: }
! 228: if (DEBUGLEVEL>=6)
! 229: fprintferr("SubCyclo: conductor:%ld\n",cnd);
! 230: avma=ltop;
! 231: return cnd;
! 232: }
! 233:
! 234: /* Calcule les orbites d'un sous-groupe de Z/nZ donne par un
! 235: * generateur ou d'un ensemble de generateur donne par un vecteur.
! 236: */
! 237: GEN
! 238: znstar_cosets(long n, long phi_n, GEN H)
! 239: {
! 240: long k;
! 241: long c = 0;
! 242: long card = group_order(H);
! 243: long index = phi_n/card;
! 244: GEN cosets = cgetg(index+1,t_VECSMALL);
! 245: gpmem_t ltop = avma;
! 246: GEN bits = bitvec_alloc(n);
! 247: for (k = 1; k <= index; k++)
! 248: {
! 249: for (c++ ; bitvec_test(bits,c) || cgcd(c,n)!=1; c++);
! 250: cosets[k]=c;
! 251: znstar_coset_bits_inplace(n, H, bits, c);
! 252: }
! 253: avma=ltop;
! 254: return cosets;
! 255: }
! 256:
! 257:
! 258: /*************************************************************************/
! 259: /** **/
! 260: /** znstar/HNF interface **/
! 261: /** **/
! 262: /*************************************************************************/
! 263:
! 264: /* Convert a true znstar output by znstar to a `small znstar'
! 265: */
! 266:
! 267: GEN
! 268: znstar_small(GEN zn)
! 269: {
! 270: GEN Z=cgetg(4,t_VEC);
! 271: Z[1]=licopy(gmael3(zn,3,1,1));
! 272: Z[2]=(long) gtovecsmall((GEN)zn[2]);
! 273: Z[3]=(long) lift((GEN)zn[3]);
! 274: return Z;
! 275: }
! 276:
! 277:
! 278: /* Compute generators for the subgroup of (Z/nZ)* given in HNF.
! 279: */
! 280: GEN
! 281: znstar_hnf_generators(GEN Z, GEN M)
! 282: {
! 283: long l = lg(M);
! 284: GEN gen=cgetg(l, t_VECSMALL);
! 285: gpmem_t ltop=avma;
! 286: GEN zgen= (GEN) Z[3];
! 287: long n = itos((GEN) Z[1]);
! 288: GEN m = stoi(n);
! 289: long j,h;
! 290: for (j = 1; j < l; j++)
! 291: {
! 292: gen[j] = 1;
! 293: for (h = 1; h < l; h++)
! 294: gen[j] = mulssmod(gen[j],
! 295: itos(powmodulo((GEN) zgen[h], gmael(M,j,h),m)),n);
! 296: }
! 297: avma=ltop;
! 298: return gen;
! 299: }
! 300:
! 301: GEN
! 302: znstar_hnf(GEN Z, GEN M)
! 303: {
! 304: return znstar_generate(itos((GEN)Z[1]),znstar_hnf_generators(Z,M));
! 305: }
! 306:
! 307: GEN
! 308: znstar_hnf_elts(GEN Z, GEN H)
! 309: {
! 310: gpmem_t ltop=avma;
! 311: GEN G=znstar_hnf(Z,H);
! 312: long n=itos((GEN)Z[1]);
! 313: GEN list=znstar_elts(n,G);
! 314: return gerepileupto(ltop,list);
! 315: }
! 316:
! 317: /*************************************************************************/
! 318: /** **/
! 319: /** subcyclo **/
! 320: /** **/
! 321: /*************************************************************************/
! 322:
! 323: static GEN gscycloconductor(GEN g, long n, long flag)
! 324: {
! 325: if (flag==2)
! 326: {
! 327: GEN V=cgetg(3,t_VEC);
! 328: V[1]=lcopy(g);
! 329: V[2]=lstoi(n);
! 330: return V;
! 331: }
! 332: return g;
! 333: }
! 334:
! 335: static long
! 336: lift_check_modulus(GEN H, long n)
! 337: {
! 338: long t=typ(H);
! 339: long h;
! 340: switch(t)
! 341: {
! 342: case t_INTMOD:
! 343: if (cmpsi(n,(GEN)H[1]))
! 344: err(talker,"wrong modulus in galoissubcyclo");
! 345: H = (GEN)H[2];
! 346: case t_INT:
! 347: h=smodis(H,n);
! 348: if (cgcd(h,n)!=1)
! 349: err(talker,"generators must be prime to conductor in galoissubcyclo");
! 350: return h;
! 351: }
! 352: err(talker,"wrong type in galoissubcyclo");
! 353: return 0;/*not reached*/
! 354: }
! 355:
! 356: GEN subcyclo_complex_bound(gpmem_t ltop, GEN V, long prec)
! 357: {
! 358: GEN pol = roots_to_pol(V,0);
! 359: GEN vec = gtovec(greal(pol));
! 360: GEN borne = ceil_safe(supnorm(vec,prec));
! 361: return gerepileupto(ltop,borne);
! 362: }
! 363:
! 364: GEN subcyclo_complex_cyclic(long n, long d, long m ,long z, long g, GEN powz, long prec)
! 365: {
! 366: GEN V=cgetg(d+1,t_VEC);
! 367: long base=1;
! 368: long i,k;
! 369: for (i=1;i<=d;i++,base=mulssmod(base,z,n))
! 370: {
! 371: gpmem_t ltop=avma;
! 372: long ex=base;
! 373: GEN s=gzero;
! 374: (void)new_chunk(2*prec + 3);
! 375: for (k=0; k<m; k++, ex = mulssmod(ex,g,n))
! 376: s=gadd(s,(GEN)powz[ex]);
! 377: avma=ltop;
! 378: V[i]=lcopy(s);
! 379: }
! 380: return V;
! 381: }
! 382:
! 383: /* Newton sums mod le. if le==NULL, works with complex instead */
! 384: GEN
! 385: subcyclo_cyclic(long n, long d, long m ,long z, long g, GEN powz, GEN le)
! 386: {
! 387: GEN V=cgetg(d+1,t_VEC);
! 388: long base=1;
! 389: long i,k;
! 390: long lle=le?lg(le)*2+1:2*lg(powz[1])+3;/*Assume dvmdii use lx+ly space*/
! 391: for (i=1;i<=d;i++,base=mulssmod(base,z,n))
! 392: {
! 393: gpmem_t ltop=avma;
! 394: long ex=base;
! 395: GEN s=gzero;
! 396: (void)new_chunk(lle); /* HACK */
! 397: for (k=0; k<m; k++, ex = mulssmod(ex,g,n))
! 398: s=gadd(s,(GEN)powz[ex]);
! 399: avma=ltop;
! 400: V[i]=le?lmodii(s,le):lcopy(s);
! 401: }
! 402: return V;
! 403: }
! 404:
! 405: struct _subcyclo_orbits_s
! 406: {
! 407: GEN powz;
! 408: GEN *s;
! 409: ulong count;
! 410: gpmem_t ltop;
! 411: };
! 412:
! 413: void
! 414: _subcyclo_orbits(struct _subcyclo_orbits_s *data, long k)
! 415: {
! 416: GEN powz = data->powz;
! 417: GEN *s = data->s;
! 418:
! 419: if (!data->count) data->ltop= avma;
! 420: *s = gadd(*s,(GEN)powz[k]);
! 421: data->count++;
! 422: if ((data->count & 0xffUL) == 0)
! 423: *s = gerepileupto(data->ltop, *s);
! 424: }
! 425:
! 426: /* Newton sums mod le. if le==NULL, works with complex instead */
! 427: GEN
! 428: subcyclo_orbits(long n, GEN H, GEN O, GEN powz, GEN le)
! 429: {
! 430: long i, d=lg(O);
! 431: GEN V=cgetg(d,t_VEC);
! 432: struct _subcyclo_orbits_s data;
! 433: long lle=le?lg(le)*2+1:2*lg(powz[1])+3;/*Assume dvmdii use lx+ly space*/
! 434: data.powz = powz;
! 435: for(i=1; i<d; i++)
! 436: {
! 437: GEN s = gzero;
! 438: gpmem_t av = avma;
! 439: (void)new_chunk(lle);
! 440: data.count = 0;
! 441: data.s = &s;
! 442: znstar_coset_func(n, H, (void (*)(void *,long)) _subcyclo_orbits,
! 443: (void *) &data, O[i]);
! 444: avma = av; /* HACK */
! 445: V[i] = le?lmodii(s,le):lcopy(s);
! 446: }
! 447: return V;
! 448: }
! 449:
! 450: GEN
! 451: subcyclo_start(long n, long d, long o, GEN borne, long *ptr_val,long *ptr_l)
! 452: {
! 453: gpmem_t av;
! 454: GEN l,le,z;
! 455: long i;
! 456: long e,val;
! 457: if (DEBUGLEVEL >= 1) (void)timer2();
! 458: l=stoi(n+1);e=1;
! 459: while(!isprime(l))
! 460: {
! 461: l=addis(l,n);
! 462: e++;
! 463: }
! 464: if (DEBUGLEVEL >= 4)
! 465: fprintferr("Subcyclo: prime l=%Z\n",l);
! 466: av=avma;
! 467: if (!borne)
! 468: {
! 469: /*Borne utilise':
! 470: Vecmax(Vec((x+o)^d)=max{binome(d,i)*o^i ;1<=i<=d}
! 471: */
! 472: i=d-(1+d)/(1+o);
! 473: borne=mulii(binome(stoi(d),i),gpowgs(stoi(o),i));
! 474: }
! 475: if (DEBUGLEVEL >= 4)
! 476: fprintferr("Subcyclo: borne=%Z\n",borne);
! 477: val=logint(shifti(borne,2),l,NULL);
! 478: avma=av;
! 479: if (DEBUGLEVEL >= 4)
! 480: fprintferr("Subcyclo: val=%ld\n",val);
! 481: le=gpowgs(l,val);
! 482: z=lift(gpowgs(gener(l),e));
! 483: z=padicsqrtnlift(gun,stoi(n),z,l,val);
! 484: if (DEBUGLEVEL >= 1)
! 485: msgtimer("padicsqrtnlift.");
! 486: *ptr_val=val;
! 487: *ptr_l=itos(l);
! 488: return gmodulcp(z,le);
! 489: }
! 490:
! 491: GEN
! 492: subcyclo_complex_roots(long n, long real, long prec)
! 493: {
! 494: GEN powz, z = exp_Ir(divrs(Pi2n(1, prec), n)); /* = e_n(1) */
! 495: long i, k = (n+3)>>1;
! 496:
! 497: powz = cgetg(n,t_VEC);
! 498: powz[1] = (long)z;
! 499: for (i=2; i<k; i++) powz[i] = lmul(z,(GEN)powz[i-1]);
! 500: if (real) /* totally real field, take real part */
! 501: {
! 502: for (i=1; i<k; i++) powz[i] = mael(powz,i,1);
! 503: for ( ; i<n; i++) powz[i] = powz[n-i];
! 504: }
! 505: else
! 506: for ( ; i<n; i++) powz[i] = lconj((GEN)powz[n-i]);
! 507: return powz;
! 508: }
! 509:
! 510: GEN
! 511: subcyclo_roots(long n, GEN zl)
! 512: {
! 513: GEN le=(GEN) zl[1];
! 514: GEN z=(GEN) zl[2];
! 515: long lle=lg(le)*3; /*Assume dvmdii use lx+ly space*/
! 516: long i;
! 517: GEN powz = cgetg(n,t_VEC);
! 518: powz[1] = (long) z;
! 519: for (i=2; i<n; i++)
! 520: {
! 521: gpmem_t av=avma;
! 522: GEN p1;
! 523: (void)new_chunk(lle); /* HACK */
! 524: p1 = mulii(z,(GEN)powz[i-1]);
! 525: avma=av;
! 526: powz[i] = lmodii(p1,le);
! 527: }
! 528: return powz;
! 529: }
! 530:
! 531: GEN
! 532: galoiscyclo(long n, long v)
! 533: {
! 534: ulong ltop=avma;
! 535: GEN grp,G;
! 536: GEN z, le;
! 537: long val,l;
! 538: GEN L;
! 539: long i,j,k;
! 540: GEN zn=znstar(stoi(n));
! 541: long card=itos((GEN) zn[1]);
! 542: GEN gen=lift((GEN)zn[3]);
! 543: GEN ord=gtovecsmall((GEN)zn[2]);
! 544: GEN elts;
! 545: z=subcyclo_start(n,card/2,2,NULL,&val,&l);
! 546: le=(GEN) z[1];
! 547: z=(GEN) z[2];
! 548: L = cgetg(1+card,t_VEC);
! 549: L[1] = (long) z;
! 550: for (j = 1, i = 1; j < lg(gen); j++)
! 551: {
! 552: int c = i * (ord[j] - 1);
! 553: for (k = 1; k <= c; k++) /* I like it */
! 554: L[++i] = (long) powmodulo((GEN)L[k],(GEN)gen[j],le);
! 555: }
! 556: G=abelian_group(ord);
! 557: elts = group_elts(G, card); /*not stack clean*/
! 558: grp = cgetg(9, t_VEC);
! 559: grp[1] = (long) cyclo(n,v);
! 560: grp[2] = lgetg(4,t_VEC);
! 561: mael(grp,2,1) = lstoi(l);
! 562: mael(grp,2,2) = lstoi(val);
! 563: mael(grp,2,3) = licopy(le);
! 564: grp[3] = lcopy(L);
! 565: grp[4] = (long) vandermondeinversemod(L, (GEN) grp[1], gun, le);
! 566: grp[5] = un;
! 567: grp[6] = lcopy(elts);
! 568: grp[7] = lcopy((GEN)G[1]);
! 569: grp[8] = lcopy((GEN)G[2]);
! 570: return gerepileupto(ltop,grp);
! 571: }
! 572:
! 573: /* Convert a bnrinit(Q,n) to a znstar(n)
! 574: * complex is set to 0 if the bnr is real and to 1 if it is complex.
! 575: * Not stack clean
! 576: */
! 577: GEN bnrtozn(GEN bnr, long *complex)
! 578: {
! 579: GEN zk;
! 580: GEN gen;
! 581: GEN cond;
! 582: long l2;
! 583: long i;
! 584: GEN p3; /* vec */
! 585: GEN res;
! 586: checkbnrgen(bnr);
! 587: zk = (GEN) bnr[5];
! 588: gen = (GEN) zk[3];
! 589: /*cond is the finite part of the conductor
! 590: * complex is the infinite part*/
! 591: cond = gcoeff(gmael3(bnr,2,1,1), 1, 1);
! 592: *complex = signe(gmael4(bnr,2,1,2,1));
! 593: l2 = lg(gen);
! 594: res= cgetg(4,t_VEC);
! 595: res[1]=zk[1];
! 596: res[2]=zk[2];
! 597: p3 = cgetg(l2, t_VEC);
! 598: for (i = 1; i < l2; ++i)
! 599: {
! 600: GEN x=(GEN) gen[i];
! 601: if (typ(x) == t_MAT)
! 602: x = gcoeff(x, 1, 1);
! 603: else if (typ(x) == t_COL)
! 604: x = (GEN) x[1];
! 605: p3[i] = (long) gmodulcp(mpabs(x), cond);
! 606: }
! 607: res[3] = (long) p3;
! 608: return res;
! 609: }
! 610:
! 611: GEN
! 612: galoissubcyclo(GEN N, GEN sg, long flag, long v)
! 613: {
! 614: gpmem_t ltop=avma,av;
! 615: GEN H, V;
! 616: long i;
! 617: GEN O;
! 618: GEN Z=NULL;
! 619: GEN B,zl,L,T,le,powz;
! 620: long val,l;
! 621: long n, cnd, complex=1;
! 622: long card, phi_n;
! 623: if (flag<0 || flag>2) err(flagerr,"galoissubcyclo");
! 624: if ( v==-1 ) v=0;
! 625: if (!sg) sg=gun;
! 626: switch(typ(N))
! 627: {
! 628: case t_INT:
! 629: n=itos(N);
! 630: if ( n<1 ) err(arither2);
! 631: break;
! 632: case t_VEC:
! 633: if (lg(N)==7)
! 634: N=bnrtozn(N,&complex);
! 635: if (lg(N)==4)
! 636: {
! 637: Z=N;
! 638: if (typ(Z[3])!=t_VEC)
! 639: err(typeer,"galoissubcyclo");
! 640: if (lg(Z[3])==1)
! 641: n=1;
! 642: else
! 643: {
! 644: if (typ(gmael(Z,3,1))!= t_INTMOD)
! 645: #ifdef NETHACK_MESSAGES
! 646: err(talker,"You have transgressed!");
! 647: #else
! 648: err(talker,"Please do not try to break PARI with ridiculously counterfeit data. Thanks!");
! 649: #endif
! 650: n=itos(gmael3(Z,3,1,1));
! 651: }
! 652: break;
! 653: }
! 654: default: /*fall through*/
! 655: err(typeer,"galoissubcyclo");
! 656: return NULL;/*Not reached*/
! 657: }
! 658: if (n==1) {avma=ltop; return polx[v];}
! 659: switch(typ(sg))
! 660: {
! 661: case t_INTMOD: case t_INT:
! 662: V=cgetg(2,t_VECSMALL);
! 663: V[1]=lift_check_modulus(sg,n);
! 664: break;
! 665: case t_VECSMALL:
! 666: V=gcopy(sg);
! 667: for (i=1;i<lg(V);i++)
! 668: if (V[i]<0)
! 669: V[i]=mulssmod(-V[i],n-1,n);
! 670: break;
! 671: case t_VEC:
! 672: case t_COL:
! 673: V=cgetg(lg(sg),t_VECSMALL);
! 674: for(i=1;i<lg(sg);i++)
! 675: V[i] = (long)lift_check_modulus((GEN)sg[i],n);
! 676: break;
! 677: case t_MAT:/*Fall through*/
! 678: {
! 679: if (lg(sg) == 1 || lg(sg) != lg(sg[1]))
! 680: err(talker,"not a HNF matrix in galoissubcyclo");
! 681: if (!Z)
! 682: err(talker,"N must be a bnrinit or a znstar if H is a matrix in galoissubcyclo");
! 683: if ( lg(Z[2]) != lg(sg) || lg(Z[3]) != lg(sg))
! 684: err(talker,"Matrix of wrong dimensions in galoissubcyclo");
! 685: V = znstar_hnf_generators(znstar_small(Z),sg);
! 686: }
! 687: break;
! 688: default:
! 689: err(typeer,"galoissubcyclo");
! 690: return NULL;/*Not reached*/
! 691: }
! 692: if (!complex) /*Add complex conjugation*/
! 693: V=vecsmall_append(V,n-1);
! 694: H = znstar_generate(n,V);
! 695: /* compute the complex/real status
! 696: * it is real iff z -> conj(z)=z^-1=z^(n-1) is in H
! 697: */
! 698: if (DEBUGLEVEL >= 6)
! 699: {
! 700: fprintferr("Subcyclo: elements:");
! 701: for (i=1;i<n;i++)
! 702: if (bitvec_test((GEN)H[3],i))
! 703: fprintferr(" %ld",i);
! 704: fprintferr("\n");
! 705: }
! 706: complex = !bitvec_test((GEN) H[3],n-1);
! 707: if (DEBUGLEVEL >= 6)
! 708: fprintferr("Subcyclo: complex=%ld\n",complex);
! 709: if (DEBUGLEVEL >= 1) (void)timer2();
! 710: cnd = znstar_conductor(n,H);
! 711: if (DEBUGLEVEL >= 1)
! 712: msgtimer("znstar_conductor");
! 713: if ( flag == 1 ) { avma=ltop; return stoi(cnd); }
! 714: if (n != cnd)
! 715: {
! 716: H = znstar_reduce_modulus(H, cnd);
! 717: n = cnd;
! 718: }
! 719: card = group_order(H);
! 720: phi_n= itos(phi(stoi(n)));
! 721: if ( card==phi_n )
! 722: {
! 723: avma=ltop;
! 724: if (flag==3) return galoiscyclo(n,v);
! 725: return gscycloconductor(cyclo(n,v),n,flag);
! 726: }
! 727: O = znstar_cosets(n, phi_n, H);
! 728: if (DEBUGLEVEL >= 1)
! 729: msgtimer("znstar_cosets");
! 730: if (DEBUGLEVEL >= 6)
! 731: fprintferr("Subcyclo: orbits=%Z\n",O);
! 732: if (DEBUGLEVEL >= 4)
! 733: fprintferr("Subcyclo: %ld orbits with %ld elements each\n",phi_n/card,card);
! 734: av=avma;
! 735: powz=subcyclo_complex_roots(n,!complex,3);
! 736: L=subcyclo_orbits(n,H,O,powz,NULL);
! 737: B=subcyclo_complex_bound(av,L,3);
! 738: zl=subcyclo_start(n,phi_n/card,card,B,&val,&l);
! 739: powz=subcyclo_roots(n,zl);
! 740: le=(GEN) zl[1];
! 741: L=subcyclo_orbits(n,H,O,powz,le);
! 742: T=FpV_roots_to_pol(L,le,v);
! 743: T=FpX_center(T,le);
! 744: return gerepileupto(ltop,gscycloconductor(T,n,flag));
! 745: }
! 746:
! 747: GEN
! 748: subcyclo(long n, long d, long v)
! 749: {
! 750: gpmem_t ltop=avma;
! 751: long o,p,al,r,g,gd;
! 752: GEN fa,G;
! 753: GEN zl,L,T,le;
! 754: long l,val;
! 755: GEN B,powz;
! 756: if (v<0) v = 0;
! 757: if (d==1) return polx[v];
! 758: if (d<=0 || n<=0) err(typeer,"subcyclo");
! 759: if ((n & 3) == 2) n >>= 1;
! 760: if (n == 1 || d >= n) err(talker,"degree does not divide phi(n) in subcyclo");
! 761: fa = decomp(stoi(n));
! 762: p = itos(gmael(fa,1,1));
! 763: al= itos(gmael(fa,2,1));
! 764: if (lg((GEN)fa[1]) > 2 || (p==2 && al>2))
! 765: err(talker,"non-cyclic case in polsubcyclo: use galoissubcyclo instead");
! 766: avma=ltop;
! 767: r = cgcd(d,n); /* = p^(v_p(d))*/
! 768: n = r*p;
! 769: o = n-r; /* = phi(n) */
! 770: if (o == d) return cyclo(n,v);
! 771: if (o % d) err(talker,"degree does not divide phi(n) in subcyclo");
! 772: o /= d;
! 773: if (p==2)
! 774: {
! 775: GEN pol = powgi(polx[v],gdeux); pol[2]=un; /* replace gzero */
! 776: return pol; /* = x^2 + 1 */
! 777: }
! 778: G=gener(stoi(n));
! 779: g=itos((GEN)G[2]);
! 780: gd=itos((GEN)gpowgs(G,d)[2]);
! 781: avma=ltop;
! 782: powz=subcyclo_complex_roots(n,(o&1)==0,3);
! 783: L=subcyclo_cyclic(n,d,o,g,gd,powz,NULL);
! 784: B=subcyclo_complex_bound(ltop,L,3);
! 785: zl=subcyclo_start(n,d,o,B,&val,&l);
! 786: le=(GEN)zl[1];
! 787: powz=subcyclo_roots(n,zl);
! 788: if (DEBUGLEVEL >= 6)
! 789: msgtimer("subcyclo_roots");
! 790: L=subcyclo_cyclic(n,d,o,g,gd,powz,le);
! 791: if (DEBUGLEVEL >= 6)
! 792: msgtimer("subcyclo_cyclic");
! 793: T=FpV_roots_to_pol(L,le,v);
! 794: if (DEBUGLEVEL >= 6)
! 795: msgtimer("roots_to_pol");
! 796: T=FpX_center(T,le);
! 797: return gerepileupto(ltop,T);
! 798: }
! 799:
! 800: GEN polsubcyclo(long n, long d, long v)
! 801: {
! 802: gpmem_t ltop=avma;
! 803: GEN L, Z=znstar(stoi(n));
! 804: /*subcyclo is twice faster but Z must be cyclic*/
! 805: if (lg(Z[2]) == 2 && divise((GEN)Z[1],stoi(d)))
! 806: {
! 807: avma=ltop;
! 808: return subcyclo(n, d, v);
! 809: }
! 810: L=subgrouplist((GEN) Z[2], _vec(stoi(d)));
! 811: if (lg(L) == 2)
! 812: return gerepileupto(ltop, galoissubcyclo(Z, (GEN) L[1], 0, v));
! 813: else
! 814: {
! 815: GEN V=cgetg(lg(L),t_VEC);
! 816: long i;
! 817: for (i=1; i< lg(V); i++)
! 818: V[i] = (long) galoissubcyclo(Z, (GEN) L[i], 0, v);
! 819: return gerepileupto(ltop,V);
! 820: }
! 821: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to subcyclo.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari-2.2 / src / basemath