File: [local] / OpenXM_contrib / pari / src / basemath / Attic / arith2.c (download)
Revision 1.1.1.1 (vendor branch), Sun Jan 9 17:35:30 2000 UTC (24 years, 6 months ago) by maekawa
Branch: PARI_GP
CVS Tags: maekawa-ipv6, VERSION_2_0_17_BETA, RELEASE_20000124, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, RELEASE_1_1_3, RELEASE_1_1_2 Changes since 1.1: +0 -0
lines
Import PARI/GP 2.0.17 beta.
|
/*********************************************************************/
/** **/
/** ARITHMETIC FUNCTIONS **/
/** (second part) **/
/** **/
/*********************************************************************/
/* $Id: arith2.c,v 1.1.1.1 1999/09/16 13:47:17 karim Exp $ */
#include "pari.h"
GEN arith_proto(long f(GEN), GEN x, int do_error);
GEN garith_proto(GEN f(GEN), GEN x, int do_error);
GEN garith_proto2gs(GEN f(GEN,long), GEN x, long y);
static long court_p[]={evaltyp(t_INT)|m_evallg(3),evalsigne(1)|evallgefint(3),0};
/***********************************************************************/
/** **/
/** PRIME NUMBERS **/
/** **/
/***********************************************************************/
GEN
prime(long n)
{
byteptr p = diffptr;
long c, prime = 0;
while (n--)
{
c = *p++; if (!c) err(primer1);
prime += c;
}
return stoi(prime);
}
GEN
primes(long n)
{
byteptr p = diffptr;
long c, prime = 0;
GEN y,z;
z = y = cgetg(n+1,t_VEC);
while (n--)
{
c = *p++; if (!c) err(primer1);
prime += c; *++z = lstoi(prime);
}
return y;
}
/* Building/Rebuilding the diffptr table. Incorporates Ilya Zakharevich's
* block sweep algorithm (see pari-dev-111, 25 April 1998). The actual work
* is done by the following two subroutines; the user entry point is the
* function initprimes() below. initprimes1() is the old algorithm, called
* when maxnum (size) is moderate.
*/
static byteptr
initprimes1(long size, long *lenp, long *lastp)
{
long k;
byteptr q,r,s,fin, p = (byteptr) gpmalloc(size+2);
memset(p, 0, size + 2); fin = p + size;
for (r=q=p,k=1; r<=fin; )
{
do { r+=k; k+=2; r+=k; } while (*++q);
for (s=r; s<=fin; s+=k) *s=1;
}
r=p; *r++=2; *r++=1; /* 2 and 3 */
for (s=q=r-1; ; s=q)
{
do q++; while (*q);
if (q>fin) break;
*r++ = (q-s) << 1;
}
*r++=0;
*lenp = r - p;
*lastp = ((s - p) << 1) + 1;
#if 0
fprintferr("initprimes1: q = %ld, len = %ld, last = %ld\n",
q, *lenp, *lastp);
#endif
return (byteptr) gprealloc(p,r-p,size + 2);
}
#define PRIME_ARENA (512 * 1024)
/* Here's the workhorse. This is recursive, although normally the first
recursive call will bottom out and invoke initprimes1() at once.
(Not static; might conceivably be useful to someone in library mode) */
byteptr
initprimes0(ulong maxnum, long *lenp, long *lastp)
{
long k, size, alloced, asize, psize, rootnum, curlow, last, remains, need;
byteptr q,r,s,fin, p, p1, fin1, plast, curdiff;
#if 0
fprintferr("initprimes0: called for maxnum = %lu\n", maxnum);
#endif
if (maxnum <= 1ul<<17) /* Arbitrary. */
return initprimes1(maxnum>>1, lenp, lastp);
maxnum |= 1; /* make it odd. */
/* Checked to be enough up to 40e6, attained at 155893 */
size = (long) (1.09 * maxnum/log((double)maxnum)) + 145;
p1 = (byteptr) gpmalloc(size);
rootnum = (long) sqrt((double)maxnum); /* cast it back to a long */
rootnum |= 1;
#if 0
fprintferr("initprimes0: rootnum = %ld\n", rootnum);
#endif
{
byteptr p2 = initprimes0(rootnum, &psize, &last); /* recursive call */
memcpy(p1, p2, psize); free(p2);
}
fin1 = p1 + psize - 1;
remains = (maxnum - rootnum) >> 1; /* number of odd numbers to check */
need = 100 * rootnum; /* Make % overhead negligeable. */
if (need < PRIME_ARENA) need = PRIME_ARENA;
if (avma - bot < need>>1) { /* need to do our own allocation */
alloced = 1; asize = need;
} else { /* scratch area is free part of PARI stack */
alloced = 0; asize = avma - bot;
}
if (asize > remains) asize = remains + 1;/* + room for a sentinel byte */
if (alloced)
p = (byteptr) gpmalloc(asize);
else
p = (byteptr) bot;
fin = p + asize - 1; /* the 0 sentinel goes at fin. */
curlow = rootnum + 2; /* We know all primes up to rootnum (odd). */
curdiff = fin1;
/* During each iteration p..fin-1 represents a range of odd
numbers. plast is a pointer which represents the last prime seen,
it may point before p..fin-1. */
plast = p - ((rootnum - last) >> 1) - 1;
while (remains)
{
if (asize > remains)
{
asize = remains + 1;
fin = p + asize - 1;
}
memset(p, 0, asize);
/* p corresponds to curlow. q runs over primediffs. */
for (q = p1+2, k = 3; q <= fin1; k += *q++)
{
/* The first odd number which is >= curlow and divisible by p
equals (if curlow > p)
the last odd number which is <= curlow + 2p - 1 and 0 (mod p)
which equals
p + the last even number which is <= curlow + p - 1 and 0 (mod p)
which equals
p + the last even number which is <= curlow + p - 2 and 0 (mod p)
which equals
p + curlow + p - 2 - (curlow + p - 2)) % 2p.
*/
long k2 = k*k - curlow;
if (k2 > 0)
{
r = p + (k2 >>= 1);
if (k2 > remains) break; /* Guard against an address wrap. */
} else
r = p - (((curlow+k-2) % (2*k)) >> 1) + k - 1;
for (s = r; s < fin; s += k) *s = 1;
}
/* now q runs over addresses corresponding to primes */
for (q = p; ; plast = q++)
{
while (*q) q++; /* use 0-sentinel at end */
if (q >= fin) break;
*curdiff++ = (q - plast) << 1;
}
plast -= asize - 1;
remains -= asize - 1;
curlow += ((asize - 1)<<1);
} /* while (remains) */
last = curlow - ((p - plast) << 1);
*curdiff++ = 0; /* sentinel */
*lenp = curdiff - p1;
*lastp = last;
if (alloced) free(p);
return (byteptr) gprealloc(p1, *lenp, size);
}
#if 0 /* not yet... GN */
/* The diffptr table will contain at least 6548 entries (up to and including
the 6547th prime, 65557, and the terminal null byte), because the isprime/
small-factor-extraction machinery wants to depend on everything up to 65539
being in the table, and we might as well go to a multiple of 4 Bytes.--GN */
void
init_tinyprimes_tridiv(byteptr p); /* in ifactor2.c */
#endif
byteptr
initprimes(long maxnum)
{
long len, last;
byteptr p;
/* The algorithm must see the next prime beyond maxnum, whence the +257. */
/* switch to unsigned: adding the 257 _could_ wrap to a negative number. */
ulong maxnum1 = ((maxnum<65302?65302:maxnum)+257ul);
if (maxnum1 > 436273000)
err(talker,"impossible to have prestored primes > 436273009");
p = initprimes0(maxnum1, &len, &last);
#if 0 /* not yet... GN */
static int build_the_tables = 1;
if (build_the_tables) { init_tinyprimes_tridiv(p); build_the_tables=0; }
#endif
return p;
}
static void
cleanprimetab(void)
{
long i,j, lp = lg(primetab);
for (i=j=1; i < lp; i++)
if (primetab[i]) primetab[j++] = primetab[i];
setlg(primetab,j);
}
/* primes is a single element or a row vector with "primes" to add to
* primetab. If the "prime" shares a non-trivial factor with another element
* in primetab, it is taken into account.
*/
GEN
addprimes(GEN prime)
{
long i,tx, av = avma, lp = lg(primetab);
GEN p1,p2;
if (!prime) return primetab;
tx = typ(prime);
if (is_vec_t(tx))
{
for (i=1; i < lg(prime); i++) addprimes((GEN) prime[i]);
return primetab;
}
if (tx != t_INT) err(typeer,"addprime");
if (is_pm1(prime)) return primetab;
if (signe(prime) < 0) prime=absi(prime);
p1=cgetg(1,t_VEC);
for (i=1; i < lp; i++)
{
p2 = mppgcd((GEN) primetab[i], prime);
if (! gcmp1(p2))
{
if (!egalii(prime,p2)) p1=concatsp(p1,p2);
p1 = concatsp(p1,divii((GEN)primetab[i],p2));
gunclone((GEN)primetab[i]); primetab[i]=0;
}
}
if (i == NUMPRTBELT+1 && lg(p1) == 1)
err(talker,"extra primetable overflows");
primetab[i] = lclone(prime); setlg(primetab, lp+1);
cleanprimetab();
if (lg(p1) > 1) addprimes(p1);
avma=av; return primetab;
}
GEN
removeprimes(GEN prime)
{
long i,tx, av = avma;
if (!prime) return primetab;
tx = typ(prime);
if (is_vec_t(tx))
{
for (i=1; i < lg(prime); i++) removeprimes((GEN) prime[i]);
return primetab;
}
if (tx != t_INT) err(typeer,"removeprimes");
if (is_pm1(prime)) return primetab;
if (signe(prime) < 0) prime=absi(prime);
for (i=1; i < lg(primetab); i++)
if (egalii((GEN) primetab[i], prime))
{
gunclone((GEN)primetab[i]); primetab[i]=0;
cleanprimetab(); avma=av; return primetab;
}
err(talker,"prime is not in primetable in removeprimes");
return NULL; /* not reached */
}
/***********************************************************************/
/** **/
/** COMPUTING THE MATRIX OF PRIME DIVISORS AND EXPONENTS **/
/** **/
/***********************************************************************/
/* Configure may/should define this to 1 if MPQS is not installed --GN */
#ifndef decomp_default_hint
#define decomp_default_hint 0
#endif
/* Some overkill removed from this (15 spsp for an integer < 2^32 ?!).
Should really revert to isprime() once the new primality tester is ready
--GN */
#define pseudoprime(p) millerrabin(p,3*lgefint(p))
/* where to stop trial dividing in factorization */
static long
tridiv_bound(GEN n, long all)
{
long size = expi(n) + 1;
if (all > 1) /* bounded factoring */
return all; /* use the given limit */
if (all == 0)
return VERYBIGINT; /* smallfact() case */
if (size <= 32)
return 16384;
else if (size <= 512)
return (size-16) << 10;
return 1L<<19; /* Rho will generally be faster above this */
}
/********** about to become obsolete --GN **********/
#if 0
/* If flag != 1, use the new MPQS code: Try to factor N using ECM if flag = 2
* and N is not too small, followed by MPQS, or just MPQS otherwise.
*/
static GEN
find_fact(GEN N, long *flag)
{
GEN f;
if (*flag == 2)
{
if ((f = pollardbrent(N))) /* assignment */
return f;
f = ellfacteur(N, 0);
if (!f)
{
/* *flag = 0; */
f = mpqs(N);
}
}
else if (*flag == 1)
{
if ((f = pollardbrent(N))) /* assignment */
return f;
f = ellfacteur(N, 1);
}
else
f = mpqs(N); /* may bail out and call ellfacteur(_,1) */
/* (this is to be changed) */
return f;
}
#endif
/***************************************************/
/* function imported from ifactor1.c */
long
ifac_decomp(GEN n, long hint);
static GEN
aux_end(GEN n, long nb)
{
GEN p,p1,p2, z = (GEN)avma;
long i;
if (n) gunclone(n);
p1=cgetg(nb+1,t_COL);
p2=cgetg(nb+1,t_COL);
for (i=nb; i; i--)
{
p2[i] = (long)z; z += lg(z);
p1[i] = (long)z; z += lg(z);
}
z[1]=(long)p1;
z[2]=(long)p2;
if (nb)
{
long av = avma;
GEN p1old = dummycopy(p1), p2old = dummycopy(p2);
p=sindexsort(p1);
for (i=1;i<=nb; i++)
{
p1[i]=p1old[p[i]];
p2[i]=p2old[p[i]];
}
avma = av;
}
return z;
}
static GEN
auxdecomp0(GEN n, long all, long hint)
{
long pp[] = { evaltyp(t_INT)|m_evallg(3), evalsigne(1)|evallgefint(3),2,0 };
long nb = 0,i,k,lim1,av,lp;
byteptr d=diffptr+1;
if (typ(n) != t_INT) err(arither1);
i=signe(n); if (!i) err(arither2);
cgetg(3,t_MAT);
if (i<0) { stoi(-1); stoi(1); nb++; }
if (is_pm1(n)) return aux_end(NULL,nb);
n = gclone(n); setsigne(n,1);
i = vali(n);
if (i)
{
stoi(2); stoi(i); nb++;
av=avma; affii(shifti(n,-i), n); avma=av;
}
if (is_pm1(n)) return aux_end(n,nb);
lim1 = tridiv_bound(n, all);
/* trial division */
while (*d && pp[2] <= lim1)
{
pp[2] += *d++;
if (mpdivis(n,pp,n))
{
nb++; k=1; while (mpdivis(n,pp,n)) k++;
icopy(pp); stoi(k);
if (is_pm1(n)) return aux_end(n,nb);
}
}
/* pp = square of biggest p tried so far */
av=avma; setlg(pp,4); affii(sqri(pp),pp); avma=av;
if (cmpii(pp,n) > 0 || pseudoprime(n))
{
nb++;
icopy(n); stoi(1); return aux_end(n,nb);
}
lp = lg(primetab); k=0;
for (i=1; i<lp; i++)
if (mpdivis(n,(GEN) primetab[i],n))
{
nb++; k=1; while (mpdivis(n,(GEN) primetab[i],n)) k++;
icopy((GEN) primetab[i]); stoi(k);
if (is_pm1(n)) return aux_end(n,nb);
}
/* test primality again, _if_ primetab made a difference */
if(k && (cmpii(pp,n) > 0 || pseudoprime(n)))
{
nb++;
icopy(n); stoi(1); return aux_end(n,nb);
}
/* now we have a large composite. Use hint as is if all==1 */
if (all!=1) hint = 15; /* turn off everything except pure powers */
nb += ifac_decomp(n, hint);
return aux_end(n, nb);
}
GEN
auxdecomp(GEN n, long all)
{
return auxdecomp0(n,all,decomp_default_hint);
}
/* see before ifac_crack() in ifactor1.c for current semantics of `hint'
(factorint's `flag') */
GEN
factorint(GEN n, long flag)
{
return auxdecomp0(n,1,flag);
}
GEN
decomp(GEN n)
{
return auxdecomp0(n,1,decomp_default_hint);
}
GEN
smallfact(GEN n)
{
return boundfact(n,0);
}
GEN
gboundfact(GEN n, long lim)
{
return garith_proto2gs(boundfact,n,lim);
}
GEN
boundfact(GEN n, long lim)
{
GEN p1,p2,p3,p4,p5,y;
long av,tetpil;
if (lim<=1) lim=0;
switch(typ(n))
{
case t_INT:
return auxdecomp(n,lim);
case t_FRACN:
av=avma; n=gred(n); /* fall through */
case t_FRAC:
if (typ(n)==t_FRAC) av=avma;
p1=auxdecomp((GEN)n[1],lim);
p2=auxdecomp((GEN)n[2],lim);
p4=concatsp((GEN)p1[1],(GEN)p2[1]);
p5=concatsp((GEN)p1[2],gneg((GEN)p2[2]));
p3=indexsort(p4);
tetpil=avma; y=cgetg(3,t_MAT);
y[1]=(long)extract(p4,p3);
y[2]=(long)extract(p5,p3);
return gerepile(av,tetpil,y);
}
err(arither1);
return NULL; /* not reached */
}
/***********************************************************************/
/** **/
/** BASIC ARITHMETIC FUNCTIONS **/
/** **/
/***********************************************************************/
/* functions imported from ifactor1.c */
long
ifac_moebius(GEN n, long hint);
long
ifac_issquarefree(GEN n, long hint);
long
ifac_omega(GEN n, long hint);
long
ifac_bigomega(GEN n, long hint);
GEN
ifac_totient(GEN n, long hint);
GEN
ifac_numdiv(GEN n, long hint);
GEN
ifac_sumdiv(GEN n, long hint);
GEN
ifac_sumdivk(GEN n, long k, long hint);
GEN
gmu(GEN n)
{
return arith_proto(mu,n,1);
}
long
mu(GEN n)
{
byteptr d = diffptr+1; /* point at 3 - 2 */
GEN p;
long av = avma, lim1, s, v;
if (typ(n) != t_INT) err(arither1);
if (!signe(n)) err(arither2);
if (is_pm1(n)) return 1;
v = vali(n);
if (v>1) return 0;
s = v ? -1 : 1;
n=absi(shifti(n,-v)); p=court_p; p[2]=2;
if (is_pm1(n)) return s;
lim1 = tridiv_bound(n,1);
while (*d && p[2] < lim1)
{
p[2] += *d++;
if (mpdivis(n,p,n))
{
if (divise(n,p)) { avma=av; return 0; }
s = -s; if (is_pm1(n)) { avma=av; return s; }
}
}
/* we normally get here with p==nextprime(PRE_TUNE), which will already
have been tested against n, so p^2 >= n (instead of p^2 > n) is
enough to guarantee n is prime */
if (cmpii(sqri(p),n) >= 0 || pseudoprime(n)) { avma=av; return -s; }
/* large composite without small factors */
v = ifac_moebius(n, decomp_default_hint);
avma=av;
return (s<0 ? -v : v); /* correct also if v==0 */
}
GEN
gissquarefree(GEN x)
{
return arith_proto(issquarefree,x,0);
}
long
issquarefree(GEN x)
{
long av=avma,lim1,tx;
GEN p;
if (gcmp0(x)) return 0;
tx=typ(x);
if (tx == t_INT)
{
long v;
byteptr d=diffptr+1;
if (is_pm1(x)) return 1;
if((v = vali(x)) > 1) return 0;
x=absi(shifti(x,-v)); p=court_p; p[2]=2;
if (is_pm1(x)) return 1;
lim1 = tridiv_bound(x,1);
while (*d && p[2] < lim1)
{
p[2] += *d++;
if (mpdivis(x,p,x))
{
if (divise(x,p)) { avma=av; return 0; }
if (is_pm1(x)) { avma=av; return 1; }
}
}
if (cmpii(sqri(p),x) >= 0 || pseudoprime(x)) { avma=av; return 1; }
v = ifac_issquarefree(x, decomp_default_hint);
avma=av;
return v;
}
if (tx != t_POL) err(typeer,"issquarefree");
p = ggcd(x, derivpol(x));
avma=av; return (lgef(p)==3);
}
GEN
gomega(GEN n)
{
return arith_proto(omega,n,1);
}
long
omega(GEN n)
{
byteptr d=diffptr+1;
GEN p;
long nb,av=avma,lim1,v;
if (typ(n) != t_INT) err(arither1);
if (!signe(n)) err(arither2);
if (is_pm1(n)) return 0;
v=vali(n);
nb = v ? 1 : 0;
n=absi(shifti(n,-v)); p=court_p; p[2]=2;
if (is_pm1(n)) return nb;
lim1 = tridiv_bound(n,1);
while (*d && p[2] < lim1)
{
p[2] += *d++;
if (mpdivis(n,p,n))
{
nb++; while (mpdivis(n,p,n)); /* empty */
if (is_pm1(n)) { avma = av; return nb; }
}
}
if (cmpii(sqri(p),n) >= 0 || pseudoprime(n)) { avma = av; return nb+1; }
/* large composite without small factors */
nb += ifac_omega(n, decomp_default_hint);
avma=av; return nb;
}
GEN
gbigomega(GEN n)
{
return arith_proto(bigomega,n,1);
}
long
bigomega(GEN n)
{
byteptr d=diffptr+1;
GEN p;
long nb,av=avma,lim1,v;
if (typ(n) != t_INT) err(arither1);
if (!signe(n)) err(arither2);
if (is_pm1(n)) return 0;
nb=v=vali(n);
n=absi(shifti(n,-v)); p=court_p; p[2]=2;
if (is_pm1(n)) { avma = av; return nb; }
lim1 = tridiv_bound(n,1);
while (*d && p[2] < lim1)
{
p[2] += *d++;
if (mpdivis(n,p,n))
{
do nb++; while (mpdivis(n,p,n));
if (is_pm1(n)) { avma = av; return nb; }
}
}
if (cmpii(sqri(p),n) >= 0 || pseudoprime(n)) { avma = av; return nb+1; }
nb += ifac_bigomega(n, decomp_default_hint);
avma=av; return nb;
}
GEN
gphi(GEN n)
{
return garith_proto(phi,n,1);
}
GEN
phi(GEN n)
{
byteptr d = diffptr+1;
GEN p,m;
long av = avma, lim1, v;
if (typ(n) != t_INT) err(arither1);
if (!signe(n)) err(arither2);
if (is_pm1(n)) return gun;
v=vali(n);
n=absi(shifti(n,-v)); p=court_p; p[2]=2;
m = (v > 1 ? shifti(gun,v-1) : stoi(1));
if (is_pm1(n)) { return gerepileupto(av,m); }
lim1 = tridiv_bound(n,1);
while (*d && p[2] < lim1)
{
court_p[2] += *d++;
if (mpdivis(n,p,n))
{
m = mulii(m,addsi(-1,p));
while (mpdivis(n,p,n)) m = mulii(m,p);
if (is_pm1(n)) { return gerepileupto(av,m); }
}
}
if (cmpii(sqri(p),n) >= 0 || pseudoprime(n))
{
m = mulii(m,addsi(-1,n));
return gerepileupto(av,m);
}
m = mulii(m, ifac_totient(n, decomp_default_hint));
return gerepileupto(av,m);
}
GEN
gnumbdiv(GEN n)
{
return garith_proto(numbdiv,n,1);
}
GEN
numbdiv(GEN n)
{
byteptr d=diffptr+1;
GEN p,m;
long l, av=avma, lim1, v;
if (typ(n) != t_INT) err(arither1);
if (!signe(n)) err(arither2);
if (is_pm1(n)) return gun;
v=vali(n);
n=absi(shifti(n,-v)); p=court_p; p[2]=2;
m = stoi(v+1);
if (is_pm1(n)) { return gerepileupto(av,m); }
lim1 = tridiv_bound(n,1);
while (*d && p[2] < lim1)
{
p[2] += *d++;
l=1; while (mpdivis(n,p,n)) l++;
m = mulsi(l,m); if (is_pm1(n)) { return gerepileupto(av,m); }
}
if(cmpii(sqri(p),n) >= 0 || pseudoprime(n))
{
m = shifti(m,1);
return gerepileupto(av,m);
}
m = mulii(m, ifac_numdiv(n, decomp_default_hint));
return gerepileupto(av,m);
}
GEN
gsumdiv(GEN n)
{
return garith_proto(sumdiv,n,1);
}
GEN
sumdiv(GEN n)
{
byteptr d=diffptr+1;
GEN p,m,m1;
long av=avma,lim1,v;
if (typ(n) != t_INT) err(arither1);
if (!signe(n)) err(arither2);
if (is_pm1(n)) return gun;
v=vali(n);
n=absi(shifti(n,-v)); p=court_p; p[2]=2;
m = (v ? addsi(-1,shifti(gun,v+1)) : stoi(1));
if (is_pm1(n)) { return gerepileupto(av,m); }
lim1 = tridiv_bound(n,1);
while (*d && p[2] < lim1)
{
p[2] += *d++;
if (mpdivis(n,p,n))
{
m1=addsi(1,p);
while (mpdivis(n,p,n)) m1=addsi(1,mulii(p,m1));
m = mulii(m1,m); if (is_pm1(n)) { return gerepileupto(av,m); }
}
}
if(cmpii(sqri(p),n) >= 0 || pseudoprime(n))
{
m = mulii(m,addsi(1,n));
return gerepileupto(av,m);
}
m = mulii(m, ifac_sumdiv(n, decomp_default_hint));
return gerepileupto(av,m);
}
GEN
gsumdivk(GEN n, long k)
{
return garith_proto2gs(sumdivk,n,k);
}
GEN
sumdivk(GEN n, long k)
{
byteptr d=diffptr+1;
GEN p,p1,m,m1,pk;
long av=avma,tetpil,k1,lim1,v;
if (!k) return numbdiv(n);
if (k==1) return sumdiv(n);
if (typ(n) != t_INT) err(arither1);
if (!signe(n)) err(arither2);
if (is_pm1(n)) return gun;
k1=k;
if (k==-1) { p1=ginv(n); m=sumdiv(n); goto fin; }
if (k<0) { k= -k; p1=gpuigs(n,k1); }
v=vali(n);
n=absi(shifti(n,-v));
p = court_p; p[2]=2; m=stoi(1);
while (v--) { m=addsi(1,shifti(m,k)); }
if (is_pm1(n)) goto fin;
lim1 = tridiv_bound(n,1);
while (*d && p[2] < lim1)
{
p[2] += *d++;
if (mpdivis(n,p,n))
{
pk=gpuigs(p,k); m1=addsi(1,pk);
while (mpdivis(n,p,n)) m1=addsi(1,mulii(pk,m1));
m = mulii(m1,m); if (is_pm1(n)) goto fin;
}
}
if(cmpii(sqri(p),n) >= 0 || pseudoprime(n))
{
pk=gpuigs(n,k);
m = mulii(m,addsi(1,pk));
goto fin;
}
m = mulii(m, ifac_sumdivk(n, k, decomp_default_hint));
fin:
if (k1>=0) return gerepileupto(av,m);
tetpil=avma; return gerepile(av,tetpil,gmul(p1,m));
}
/***********************************************************************/
/** **/
/** MISCELLANEOUS ARITHMETIC FUNCTIONS **/
/** (all of these ultimately depend on auxdecomp()) **/
/** **/
/***********************************************************************/
GEN
divisors(GEN n)
{
long tetpil,av=avma,i,j,l;
GEN *d,*t,*t1,*t2,*t3, nbdiv,e;
if (typ(n) != t_MAT || lg(n) != 3) n = auxdecomp(n,1);
e = (GEN) n[2], n = (GEN) n[1]; l = lg(n);
if (l>1 && signe(n[1]) < 0) { e++; n++; l--; } /* skip -1 */
nbdiv = gun;
for (i=1; i<l; i++)
{
e[i] = itos((GEN)e[i]);
nbdiv = mulis(nbdiv,1+e[i]);
}
d = t = (GEN*) cgetg(itos(nbdiv)+1,t_VEC);
*++d = gun;
for (i=1; i<l; i++)
for (t1=t,j=e[i]; j; j--,t1=t2)
for (t2=d,t3=t1; t3<t2; )
*++d = mulii(*++t3, (GEN)n[i]);
tetpil=avma; return gerepile(av,tetpil,sort((GEN)t));
}
GEN
core(GEN n)
{
long av=avma,tetpil,i;
GEN fa,p1,p2,res=gun;
fa=auxdecomp(n,1); p1=(GEN)fa[1]; p2=(GEN)fa[2];
for (i=1; i<lg(p1); i++)
if (mod2((GEN)p2[i])) res = mulii(res,(GEN)p1[i]);
tetpil=avma; return gerepile(av,tetpil,icopy(res));
}
GEN
core2(GEN n)
{
long av=avma,tetpil,i;
GEN fa,p1,p2,p3,res=gun,res2=gun,y;
fa=auxdecomp(n,1); p1=(GEN)fa[1]; p2=(GEN)fa[2];
for (i=1; i<lg(p1); i++)
{
p3=(GEN)p2[i];
if (mod2(p3)) res=mulii(res,(GEN)p1[i]);
if (!gcmp1(p3)) res2=mulii(res2,gpui((GEN)p1[i],shifti(p3,-1),0));
}
tetpil=avma; y=cgetg(3,t_VEC);
y[1]=(long)icopy(res); y[2]=(long)icopy(res2);
return gerepile(av,tetpil,y);
}
GEN
core0(GEN n,long flag)
{
return flag? core2(n): core(n);
}
GEN
coredisc(GEN n)
{
long av=avma,tetpil,r;
GEN p1=core(n);
r=mod4(p1); if (signe(p1)<0) r = 4-r;
if (r==1 || r==4) return p1;
tetpil=avma; return gerepile(av,tetpil,shifti(p1,2));
}
GEN
coredisc2(GEN n)
{
long av=avma,tetpil,r;
GEN y,p1,p2=core2(n);
p1=(GEN)p2[1]; r=mod4(p1); if (signe(p1)<0) r=4-r;
if (r==1 || r==4) return p2;
tetpil=avma; y=cgetg(3,t_VEC);
y[1]=lshifti(p1,2); y[2]=lmul2n((GEN)p2[2],-1);
return gerepile(av,tetpil,y);
}
GEN
coredisc0(GEN n,long flag)
{
return flag? coredisc2(n): coredisc(n);
}
/***********************************************************************/
/** **/
/** BINARY QUADRATIC FORMS **/
/** **/
/***********************************************************************/
GEN
qf_disc(GEN x, GEN y, GEN z)
{
if (!y) { y=(GEN)x[2]; z=(GEN)x[3]; x=(GEN)x[1]; }
return subii(sqri(y), shifti(mulii(x,z),2));
}
static GEN
qf_create(GEN x, GEN y, GEN z, long s)
{
GEN t;
if (typ(x)!=t_INT || typ(y)!=t_INT || typ(z)!=t_INT) err(typeer,"Qfb");
if (!s)
{
long av=avma; s = signe(qf_disc(x,y,z)); avma=av;
if (!s) err(talker,"zero discriminant in Qfb");
}
t = (s>0)? cgetg(5,t_QFR): cgetg(4,t_QFI);
t[1]=(long)icopy(x); t[2]=(long)icopy(y); t[3]=(long)icopy(z);
return t;
}
GEN
qfi(GEN x, GEN y, GEN z)
{
return qf_create(x,y,z,-1);
}
GEN
qfr(GEN x, GEN y, GEN z, GEN d)
{
GEN t = qf_create(x,y,z,1);
if (typ(d) != t_REAL)
err(talker,"Shanks distance should be a t_REAL (in qfr)");
t[4]=lrcopy(d); return t;
}
GEN
Qfb0(GEN x, GEN y, GEN z, GEN d, long prec)
{
GEN t = qf_create(x,y,z,0);
if (lg(t)==4) return t;
if (!d) d = gzero;
if (typ(d) == t_REAL)
t[4] = lrcopy(d);
else
{ t[4]=lgetr(prec); gaffect(d,(GEN)t[4]); }
return t;
}
/* composition */
static void
comp_gen(GEN z,GEN x,GEN y)
{
GEN s,n,d,d1,x1,x2,y1,y2,v1,v2,v3,b3,c3,m,p1,r;
s=shifti(addii((GEN)x[2],(GEN)y[2]),-1);
n=subii((GEN)y[2],s);
d = bezout((GEN)y[1],(GEN)x[1],&y1,&x1);
d1 = bezout(s,d,&x2,&y2);
if (gcmp1(d1))
{
v3 = (GEN)x[1];
v2 = (GEN)y[1];
}
else
{
v1=divii((GEN)x[1],d1);
v2=divii((GEN)y[1],d1);
v3 = mulii(v1, mppgcd(d1,mppgcd((GEN)x[3],mppgcd((GEN)y[3],n))));
}
m = addii(mulii(mulii(y1,y2),n), mulii((GEN)y[3],x2));
setsigne(m,-signe(m));
r=modii(m,v3); p1=mulii(v2,r); b3=shifti(p1,1);
c3=addii(mulii((GEN)y[3],d1),mulii(r,addii((GEN)y[2],p1)));
z[1]=lmulii(v3,v2);
z[2]=laddii((GEN)y[2],b3);
z[3]=ldivii(c3,v3);
}
static GEN
compimag0(GEN x, GEN y, int raw)
{
long tx = typ(x), av = avma, tetpil;
GEN z;
if (typ(y) != tx || tx!=t_QFI) err(typeer,"composition");
if (cmpii((GEN)x[1],(GEN)y[1]) > 0) { z=x; x=y; y=z; }
z=cgetg(4,t_QFI); comp_gen(z,x,y); tetpil=avma;
return gerepile(av,tetpil,raw? gcopy(z): redimag(z));
}
static GEN
compreal0(GEN x, GEN y, int raw)
{
long tx = typ(x), av = avma, tetpil;
GEN z;
if (typ(y) != tx || tx!=t_QFR) err(typeer,"composition");
z=cgetg(5,t_QFR); comp_gen(z,x,y);
z[4]=laddrr((GEN)x[4],(GEN)y[4]); tetpil=avma;
return gerepile(av,tetpil, raw? gcopy(z): redreal(z));
}
GEN
compreal(GEN x, GEN y) { return compreal0(x,y,0); }
GEN
comprealraw(GEN x, GEN y) { return compreal0(x,y,1); }
GEN
compimag(GEN x, GEN y) { return compimag0(x,y,0); }
GEN
compimagraw(GEN x, GEN y) { return compimag0(x,y,1); }
GEN
compraw(GEN x, GEN y)
{
return (typ(x)==t_QFI)? compimagraw(x,y): comprealraw(x,y);
}
static void
sq_gen(GEN z, GEN x)
{
GEN d1,x2,y2,v1,v3,b3,c3,m,p1,r;
d1 = bezout((GEN)x[2],(GEN)x[1],&x2,&y2);
if (gcmp1(d1)) v1 = v3 = (GEN)x[1];
else
{
v1 = divii((GEN)x[1],d1);
v3 = mulii(v1,mppgcd(d1,(GEN)x[3]));
}
m=mulii((GEN)x[3],x2); setsigne(m,-signe(m));
r=modii(m,v3); p1=mulii(v1,r); b3=shifti(p1,1);
c3=addii(mulii((GEN)x[3],d1),mulii(r,addii((GEN)x[2],p1)));
z[1]=lmulii(v3,v1);
z[2]=laddii((GEN)x[2],b3);
z[3]=ldivii(c3,v3);
}
GEN
sqcompimag0(GEN x, long raw)
{
long av = avma, tetpil;
GEN z = cgetg(4,t_QFI);
if (typ(x)!=t_QFI) err(typeer,"composition");
sq_gen(z,x); tetpil=avma;
return gerepile(av,tetpil,raw? gcopy(z): redimag(z));
}
static GEN
sqcompreal0(GEN x, long raw)
{
long av = avma, tetpil;
GEN z = cgetg(5,t_QFR);
if (typ(x)!=t_QFR) err(typeer,"composition");
sq_gen(z,x); z[4]=lshiftr((GEN)x[4],1); tetpil=avma;
return gerepile(av,tetpil,raw? gcopy(z): redreal(z));
}
GEN
sqcompreal(GEN x) { return sqcompreal0(x,0); }
GEN
sqcomprealraw(GEN x) { return sqcompreal0(x,1); }
GEN
sqcompimag(GEN x) { return sqcompimag0(x,0); }
GEN
sqcompimagraw(GEN x) { return sqcompimag0(x,1); }
GEN
real_unit_form(GEN x)
{
long av = avma, tetpil,prec;
GEN p1,p2, y = cgetg(5,t_QFR);
y[1]=un; p1=qf_disc(x,NULL,NULL); p2=racine(p1);
/* we know that p1 and p2 are non-zero */
if (mod2(p1) != mod2(p2)) p2=addsi(-1,p2);
y[2]=(long)p2;
y[3]=lshifti(subii(sqri(p2),p1),-2);
prec = precision((GEN)x[4]);
if (!prec)
err(talker,"not a type real in 4th component of a t_QFR");
y[4]=(long)realzero(prec); tetpil=avma;
return gerepile(av,tetpil,gcopy(y));
}
GEN
imag_unit_form(GEN x)
{
long av,tetpil;
GEN p1,p2, y = cgetg(4,t_QFI);
y[1]=un;
y[2]=mpodd((GEN)x[2])? un: zero;
av=avma; p1=mulii((GEN)x[1],(GEN)x[3]);
p2=shifti(sqri((GEN)x[2]),-2); tetpil=avma;
y[3]=lpile(av,tetpil,subii(p1,p2));
return y;
}
GEN
powrealraw(GEN x, long n)
{
long av = avma, m;
GEN y;
if (typ(x) != t_QFR)
err(talker,"not a real quadratic form in powreal");
if (!n) return real_unit_form(x);
if (n== 1) return gcopy(x);
if (n==-1) return ginv(x);
y = NULL; m=labs(n);
for (; m>1; m>>=1)
{
if (m&1) y = y? comprealraw(y,x): x;
x=sqcomprealraw(x);
}
y = y? comprealraw(y,x): x;
if (n<0) y = ginv(y);
return gerepileupto(av,y);
}
GEN
powimagraw(GEN x, long n)
{
long av = avma, m;
GEN y;
if (typ(x) != t_QFI)
err(talker,"not an imaginary quadratic form in powimag");
if (!n) return imag_unit_form(x);
if (n== 1) return gcopy(x);
if (n==-1) return ginv(x);
y = NULL; m=labs(n);
for (; m>1; m>>=1)
{
if (m&1) y = y? compimagraw(y,x): x;
x=sqcompimagraw(x);
}
y = y? compimagraw(y,x): x;
if (n<0) y = ginv(y);
return gerepileupto(av,y);
}
GEN
powraw(GEN x, long n)
{
return (typ(x)==t_QFI)? powimagraw(x,n): powrealraw(x,n);
}
/* composition: Shanks' NUCOMP & NUDUPL */
/* l = floor((|d|/4)^(1/4)) */
GEN
nucomp(GEN x, GEN y, GEN l)
{
long av=avma,tetpil,cz;
GEN a,a1,a2,b,b2,d,d1,e,g,n,p1,p2,p3,q1,q2,q3,q4,s,t2,t3,u,u1,v,v1,v2,v3,z;
if (x==y) return nudupl(x,l);
if (typ(x) != t_QFI || typ(y) != t_QFI)
err(talker,"not an imaginary quadratic form in nucomp");
if (cmpii((GEN)x[1],(GEN)y[1]) < 0) { s=x; x=y; y=s; }
s=shifti(addii((GEN)x[2],(GEN)y[2]),-1); n=subii((GEN)y[2],s);
a1=(GEN)x[1]; a2=(GEN)y[1]; d=bezout(a2,a1,&u,&v);
if (gcmp1(d)) { a=negi(gmul(u,n)); d1=d; }
else
if (divise(s,d))
{
a=negi(gmul(u,n)); d1=d; a1=divii(a1,d1);
a2=divii(a2,d1); s=divii(s,d1);
}
else
{
d1=bezout(s,d,&u1,&v1);
if (!gcmp1(d1))
{
a1=divii(a1,d1); a2=divii(a2,d1);
s=divii(s,d1); d=divii(d,d1);
}
p1=resii((GEN)x[3],d); p2=resii((GEN)y[3],d);
p3=modii(negi(mulii(u1,addii(mulii(u,p1),mulii(v,p2)))),d);
a=subii(mulii(p3,divii(a1,d)),mulii(u,divii(n,d)));
}
a=modii(a,a1); p1=subii(a1,a); if (cmpii(a,p1)>0) a=negi(p1);
v=gzero; d=a1; v2=gun; v3=a;
for (cz=0; absi_cmp(v3,l) > 0; cz++)
{
p1=dvmdii(d,v3,&t3); t2=subii(v,mulii(p1,v2));
v=v2; d=v3; v2=t2; v3=t3;
}
z=cgetg(4,t_QFI);
if (!cz)
{
q1=mulii(a2,v3); q2=addii(q1,n);
g=divii(addii(mulii(v3,s),(GEN)y[3]),d);
z[1]=lmulii(d,a2);
z[2]=laddii(shifti(q1,1),(GEN)y[2]);
z[3]=laddii(mulii(v3,divii(q2,d)), mulii(g,d1));
}
else
{
if (cz&1) { v3=negi(v3); v2=negi(v2); }
b=divii(addii(mulii(a2,d),mulii(n,v)),a1);
q1=mulii(b,v3); q2=addii(q1,n);
e=divii(addii(mulii(s,d),mulii((GEN)y[3],v)),a1);
q3=mulii(e,v2); q4=subii(q3,s);
g=divii(q4,v); b2=addii(q3,q4);
if (!gcmp1(d1)) { v2=mulii(d1,v2); v=mulii(d1,v); b2=mulii(d1,b2); }
z[1]=laddii(mulii(d,b),mulii(e,v));
z[2]=laddii(b2,addii(q1,q2));
z[3]=laddii(mulii(v3,divii(q2,d)), mulii(g,v2));
}
tetpil=avma; return gerepile(av,tetpil,redimag(z));
}
GEN
nudupl(GEN x, GEN l)
{
long av=avma,tetpil,cz;
GEN u,v,d,d1,p1,a,b,c,b2,z,v2,v3,t2,t3,e,g;
if (typ(x) != t_QFI)
err(talker,"not an imaginary quadratic form in nudupl");
d1=bezout((GEN)x[2],(GEN)x[1],&u,&v);
a=divii((GEN)x[1],d1); b=divii((GEN)x[2],d1);
c=modii(negi(mulii(u,(GEN)x[3])),a); p1=subii(a,c);
if (cmpii(c,p1)>0) c=negi(p1);
v=gzero; d=a; v2=gun; v3=c;
for (cz=0; absi_cmp(v3,l) > 0; cz++)
{
p1=dvmdii(d,v3,&t3); t2=subii(v,mulii(p1,v2));
v=v2; d=v3; v2=t2; v3=t3;
}
z=cgetg(4,t_QFI);
if (!cz)
{
g=divii(addii(mulii(v3,b),(GEN)x[3]),d);
z[1]=(long)sqri(d);
z[2]=laddii((GEN)x[2],shifti(mulii(d,v3),1));
z[3]=laddii(sqri(v3),mulii(g,d1));
}
else
{
if (cz&1) { v=negi(v); d=negi(d); }
e=divii(addii(mulii((GEN)x[3],v),mulii(b,d)),a);
g=divii(subii(mulii(e,v2),b),v);
b2=addii(mulii(e,v2),mulii(v,g));
if (!gcmp1(d1)) { v2=mulii(d1,v2); v=mulii(d1,v); b2=mulii(d1,b2); }
z[1]=laddii(sqri(d),mulii(e,v));
z[2]=laddii(b2,shifti(mulii(d,v3),1));
z[3]=laddii(sqri(v3),mulii(g,v2));
}
tetpil=avma; return gerepile(av,tetpil,redimag(z));
}
GEN
nupow(GEN x, GEN n)
{
long av,tetpil,i,j;
unsigned long m;
GEN y,l;
if (typ(n) != t_INT) err(talker,"not an integer exponent in nupow");
if (gcmp1(n)) return gcopy(x);
av=avma; y = imag_unit_form(x);
if (!signe(n)) return y;
l = racine(shifti(racine((GEN)y[3]),1));
for (i=lgefint(n)-1; i>2; i--)
for (m=n[i],j=0; j<BITS_IN_LONG; j++,m>>=1)
{
if (m&1) y=nucomp(y,x,l);
x=nudupl(x,l);
}
for (m=n[2]; m>1; m>>=1)
{
if (m&1) y=nucomp(y,x,l);
x=nudupl(x,l);
}
tetpil=avma; y=nucomp(y,x,l);
if (signe(n)<0 && !egalii((GEN)y[1],(GEN)y[2])
&& !egalii((GEN)y[1],(GEN)y[3]))
setsigne(y[2],-signe(y[2]));
return gerepile(av,tetpil,y);
}
/* reduction */
static GEN
abs_dvmdii(GEN b, GEN p1, GEN *pt, long s)
{
if (s<0) setsigne(b, 1); p1 = dvmdii(b,p1,pt);
if (s<0) setsigne(b,-1); return p1;
}
static GEN
rhoimag0(GEN x, long *flag)
{
GEN p1,b,d,z;
long fl, s = signe(x[2]);
fl = cmpii((GEN)x[1], (GEN)x[3]);
if (fl <= 0)
{
long fg = absi_cmp((GEN)x[1], (GEN)x[2]);
if (fg >= 0)
{
*flag = (s<0 && (!fl || !fg))? 2 /* set x[2] = negi(x[2]) in caller */
: 1;
return x;
}
}
p1=shifti((GEN)x[3],1); d = abs_dvmdii((GEN)x[2],p1,&b,s);
if (s>=0)
{
setsigne(d,-signe(d));
if (cmpii(b,(GEN)x[3])<=0) setsigne(b,-signe(b));
else { d=addsi(-1,d); b=subii(p1,b); }
}
else if (cmpii(b,(GEN)x[3])>=0) { d=addsi(1,d); b=subii(b,p1); }
z=cgetg(4,t_QFI);
z[1] = x[3];
z[3] = laddii((GEN)x[1], mulii(d,shifti(subii((GEN)x[2],b),-1)));
if (signe(b)<0 && !fl) setsigne(b,1);
z[2] = (long)b; *flag = 0; return z;
}
/* if sqrtD != NULL, compute Shanks' distance as well */
static GEN
rhoreal0(GEN x, GEN D, GEN isqrtD, GEN sqrtD)
{
long s = signe(x[3]);
GEN p1,p2, z = cgetg(5,t_QFR);
if (!sqrtD)
z[4] = x[4];
else
{
p1 = divrr(addir((GEN)x[2],sqrtD), subir((GEN)x[2],sqrtD));
if (signe(p1)<0) setsigne(p1,1); /* p1 = |(b+sqrtD) / (b-sqrtD)| */
z[4] = laddrr(shiftr(mplog(p1),-1), (GEN) x[4]);
}
z[1] = x[3]; setsigne(z[1],1);
p1=shifti((GEN)z[1],1);
p2 = (cmpii(isqrtD,(GEN)z[1]) >= 0)? isqrtD: (GEN)z[1];
p2 = divii(addii(p2,(GEN)x[2]),p1);
z[2] = lsubii(mulii(p2,p1), (GEN)x[2]);
setsigne(z[1],s);
p1=shifti(subii(sqri((GEN)z[2]),D),-2);
z[3] = ldivii(p1,(GEN)z[1]); return z;
}
static GEN
redreal0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD)
{
long av=avma,l,step, use_d;
if (typ(x) != t_QFR) err(talker,"not a real quadratic form in redreal");
switch(flag)
{
case 0: use_d=1; step=0; break;
case 1: use_d=1; step=1; break;
case 2: use_d=0; step=0; break;
case 3: use_d=0; step=1; break;
default: err(flagerr,"qfbred");
}
if (!D) D = qf_disc(x,NULL,NULL);
else if (typ(D) != t_INT) err(arither1);
if (!isqrtD || (use_d && !sqrtD))
{
l=max(lg(x[4]),3);
l=max(l,((BITS_IN_LONG-1-expo(x[4]))>>TWOPOTBITS_IN_LONG)+2);
sqrtD=gsqrt(D,l); isqrtD=mptrunc(sqrtD);
}
else
{
if (typ(isqrtD) != t_INT) err(arither1);
if (sqrtD)
{
l=typ(sqrtD); if (!is_intreal_t(l)) err(arither1);
}
}
if (step)
x = rhoreal0(x,D,isqrtD,sqrtD);
else
for(;;)
{
if (signe(x[2]) > 0 && cmpii((GEN)x[2],isqrtD) <= 0 )
{
GEN p1=subii(isqrtD, shifti(absi((GEN)x[1]),1));
l = absi_cmp((GEN)x[2],p1);
if (l>0 || (l==0 && signe(p1)<0)) break;
}
x = rhoreal0(x,D,isqrtD,sqrtD);
}
l=avma; return gerepile(av,l,gcopy(x));
}
GEN
redimag(GEN x)
{
long av=avma, tetpil, fl;
do x = rhoimag0(x, &fl); while (fl == 0);
tetpil=avma; x = gerepile(av,tetpil,gcopy(x));
if (fl == 2) setsigne(x[2], -signe(x[2]));
return x;
}
GEN
redreal(GEN x)
{
return redreal0(x,0,NULL,NULL,NULL);
}
GEN
rhoreal(GEN x)
{
return redreal0(x,1,NULL,NULL,NULL);
}
GEN
redrealnod(GEN x, GEN isqrtD)
{
return redreal0(x,2,NULL,isqrtD,NULL);
}
GEN
rhorealnod(GEN x, GEN isqrtD)
{
return redreal0(x,3,NULL,isqrtD,NULL);
}
GEN
qfbred0(GEN x, long flag, GEN D, GEN isqrtD, GEN sqrtD)
{
long tx=typ(x),av,tetpil,fl;
if (!is_qf_t(tx)) err(talker,"not a quadratic form in qfbred");
if (!D) D = qf_disc(x,NULL,NULL);
switch(signe(D))
{
case 1 :
return redreal0(x,flag,D,isqrtD,sqrtD);
case -1:
if (!flag) return redimag(x);
if (flag !=1) err(flagerr,"qfbred");
av=avma; x = rhoimag0(x,&fl); tetpil=avma;
x = gerepile(av,tetpil,gcopy(x));
if (fl == 2) setsigne(x[2], -signe(x[2]));
return x;
}
err(redpoler,"qfbred");
return NULL; /* not reached */
}
GEN
primeform(GEN x, GEN p, long prec)
{
long av,tetpil,s=signe(x);
GEN y,b;
if (typ(x) != t_INT || !s) err(arither1);
if (typ(p) != t_INT || signe(p) <= 0) err(arither1);
if (s < 0)
{
y = cgetg(4,t_QFI); s = 8-mod8(x);
}
else
{
y = cgetg(5,t_QFR); s = mod8(x);
y[4]=(long)realzero(prec);
}
switch(s&3)
{
case 2: case 3: err(funder2,"primeform");
}
y[1] = (long)icopy(p); av=avma;
if (egalii(p,gdeux))
{
switch(s)
{
case 0: y[2]=zero; break;
case 8: s=0; y[2]=zero; break;
case 1: s=1; y[2]=un; break;
case 4: s=4; y[2]=deux; break;
default: err(sqrter5);
}
b=subsi(s,x); tetpil=avma; b=shifti(b,-3);
}
else
{
b = mpsqrtmod(x,p); if (!b) err(sqrter5);
if (mod2(b) == mod2(x)) y[2] = (long)b;
else
{ tetpil = avma; y[2] = lpile(av,tetpil,subii(p,b)); }
av=avma; b=shifti(subii(sqri((GEN)y[2]),x),-2);
tetpil=avma; b=divii(b,p);
}
y[3]=lpile(av,tetpil,b); return y;
}