PLASMA  2.4.5
PLASMA - Parallel Linear Algebra for Scalable Multi-core Architectures
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
Macros | Functions
gkkleader.c File Reference
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "common.h"
#include "primes.h"
#include "gkkleader.h"
Include dependency graph for gkkleader.c:

Go to the source code of this file.

Macros

#define GCAND_SIZE   46

Functions

int GKK_doublingtable (int x, int m, int emax, int *dt)
int GKK_modpow (int *dt, int e, int m)
int GKK_primroot (int p, int e, primedec_t *pr_p1, int t_p1)
int GKK_multorder (int n, int p, int e, int pe, primedec_t *pr_p1, int t_p1)
void GKK_prepare (int q, int n, primedec_t *pr_q, int *t, primedec_t **pr_pi1, int *t_pi1, int *gi, int *Nif, int *Kif, int *dif)
void GKK_L (int t, primedec_t *pr_q, int *fi, int *Kif, int *dif, int *Li, int *diLi, int *cl, int *nl)
void GKK_precompute_terms (int q, primedec_t *pr_q, int t, int *gi, int *diLi, int *rp, int *Mg, int nMg)
void GKK_BalanceLoad (int thrdnbr, int *Tp, int *leaders, int nleaders, int L)
void GKK_output_tables (int m, int n, int q, primedec_t *pr_q, int t, int *gi, int *Nif, int *Kif, int *dif)
int GKK_getLeaderNbr (int me, int ne, int *nleaders, int **leaders)

Detailed Description

PLASMA InPlaceTransformation module PLASMA is a software package provided by Univ. of Tennessee, Univ. of California Berkeley and Univ. of Colorado Denver

This work is the implementation of an inplace transformation based on the GKK algorithm by Gustavson, Karlsson, Kagstrom and its fortran implementation.

Version:
2.4.5
Author:
Mathieu Faverge
Date:
2010-11-15

Definition in file gkkleader.c.

Macro Definition Documentation

#define GCAND_SIZE   46

GKK_primroot finds a primitive root of p^e, given the prime decomposition of p-1.

See also:
GKK_doublingtable
GKK_modpow
Parameters:
[in]p
[in]e
[in]pr_p1Prime decomposition of p - 1
[in]t_p1Number of primes in prime decomposition of p-1
Returns:
Return values:
Aprimitive root of p^e

Definition at line 149 of file gkkleader.c.

Function Documentation

void GKK_BalanceLoad ( int  thrdnbr,
int *  Tp,
int *  leaders,
int  nleaders,
int  L 
)

GKK_BalanceLoad balances the load by splitting across some cycles.

Parameters:
[in]thrdnbrnumber of thread working on cycles
[in,out]TpArray of size thrdnbr. Work load for each thread.
[in,out]leadersArray of size nleaders containing leaders, cycles length and owners
[in]nleadersNumber of leaders stored multiplied by 3
[in]LChunk size to compute the load balancing

Definition at line 636 of file gkkleader.c.

References IMBALANCE_THRESHOLD, L, maxval(), min, and sum().

{
int i, j, n1, nel;
double sumTi, maxTi;
/* quick return */
if( thrdnbr == 1 )
return;
sumTi = (double)sum(thrdnbr, Tp);
maxTi = (double)maxval(thrdnbr, Tp);
/* check if the load is too imbalanced */
if( (1.0 - (sumTi / ((double)thrdnbr * maxTi))) > IMBALANCE_THRESHOLD ) {
/* split across cycle */
for (i=0; i<nleaders; i+=3) {
j = leaders[i+2];
nel = leaders[i+1];
if( (nel >= thrdnbr) &&
((double)Tp[j] > (sumTi / ((double)thrdnbr * (1.0 - IMBALANCE_THRESHOLD)))) ) {
/* split */
n1 = (nel + thrdnbr - 1) / thrdnbr;
Tp[j] = Tp[j] - nel * L;
for (j=0; j<thrdnbr; j++)
Tp[j] = Tp[j] + min(n1, nel - n1*(j-1));
maxTi = (double)maxval(thrdnbr, Tp);
leaders[i+2] = -2;
}
}
}
}

Here is the call graph for this function:

Here is the caller graph for this function:

int GKK_doublingtable ( int  x,
int  m,
int  emax,
int *  dt 
)

GKK_doublingtable Computes a table of x^i mod m for i = 1, 2, 4, 8, ..., emax

Parameters:
[in]x
[in]m
[in]emaxMaximum power to compute
[out]dtOn exit, dt contains (x^1 mod m, x^2 mod m, x^4 mod m, ..., x^emax mod m)

Definition at line 47 of file gkkleader.c.

References PLASMA_ERR_OUT_OF_RESOURCES, plasma_error(), PLASMA_SUCCESS, and PWR_MAXSIZE.

{
int64_t y8, m8;
int i, numbits;
int z;
/* Convert in int64 */
m8 = m;
/* find the number of bits in emax */
numbits = 0;
z = emax;
while( z > 0 ) {
numbits = numbits + 1;
z = z / 2;
}
if (numbits > PWR_MAXSIZE) {
plasma_error(__func__, "PWR_MAXSIZE too small");
return PLASMA_ERR_OUT_OF_RESOURCES;
}
/* compute doubling table */
y8 = x;
dt[0] = x;
for (i=1; i < numbits; i++){
y8 = (y8*y8) % m8;
dt[i] = (int)y8;
}
return PLASMA_SUCCESS;
}

Here is the call graph for this function:

Here is the caller graph for this function:

int GKK_getLeaderNbr ( int  me,
int  ne,
int *  nleaders,
int **  leaders 
)

Definition at line 753 of file gkkleader.c.

References primedec::e, GKK_L(), GKK_output_tables(), GKK_precompute_terms(), GKK_prepare(), primedec::p, primedec::pe, PLASMA_SUCCESS, PRIME_MAXSIZE, PWR_MAXSIZE, and SIZE_MG.

{
primedec_t pr_q[PRIME_MAXSIZE];
primedec_t **pr_pi1 = NULL;
int *ptr;
int *t_pi1 = NULL; /* Number of primes in prime decomposition of each pi */
int *gi = NULL; /* Prime roots of pi^ei */
int *fi = NULL; /* Copy of ei-1 because it's the indice in Kif */
int *pifi = NULL; /* Copy of pi^ei */
int *rp = NULL; /* */
int *Li = NULL; /* gcd(Ki, lcm(Kh .. Kt)) with h=0/1 (q even or odd) */
int *si = NULL; /* */
int *diLi = NULL; /* di * Li */
int *Mg = NULL; /* */
int *vMg = NULL; /* */
int *Kif = NULL; /* */
int *Nif = NULL; /* */
int *dif = NULL; /* */
int i, j, s;
int q, t = 0, ndiv, ht, div;
int cl, nl, nlead;
ptr = (int *)malloc( ((3*PWR_MAXSIZE + 8)*PRIME_MAXSIZE + 4 + 2*SIZE_MG) * sizeof(int) );
t_pi1 = ptr; /* plasma_malloc(t_pi1, PRIME_MAXSIZE, int); */
gi = t_pi1 + PRIME_MAXSIZE; /* plasma_malloc(gi , PRIME_MAXSIZE, int); */
fi = gi + PRIME_MAXSIZE; /* plasma_malloc(fi , PRIME_MAXSIZE, int); */
pifi = fi + PRIME_MAXSIZE; /* plasma_malloc(pifi , PRIME_MAXSIZE, int); */
rp = pifi + PRIME_MAXSIZE; /* plasma_malloc(rp , PRIME_MAXSIZE+1, int); */
Li = rp + PRIME_MAXSIZE+1; /* plasma_malloc(Li , PRIME_MAXSIZE+1, int); */
si = Li + PRIME_MAXSIZE+1; /* plasma_malloc(si , PRIME_MAXSIZE+1, int); */
diLi = si + PRIME_MAXSIZE+1; /* plasma_malloc(diLi , PRIME_MAXSIZE+1, int); */
Kif = diLi + PRIME_MAXSIZE+1; /* plasma_malloc(Kif , PRIME_MAXSIZE*PWR_MAXSIZE, int); */
Nif = Kif + PRIME_MAXSIZE*PWR_MAXSIZE; /* plasma_malloc(Nif , PRIME_MAXSIZE*PWR_MAXSIZE, int); */
dif = Nif + PRIME_MAXSIZE*PWR_MAXSIZE; /* plasma_malloc(dif , PRIME_MAXSIZE*PWR_MAXSIZE, int); */
Mg = dif + PRIME_MAXSIZE*PWR_MAXSIZE; /* plasma_malloc(Mg , SIZE_MG, int); */
vMg = Mg + SIZE_MG; /* plasma_malloc(vMg , SIZE_MG, int); */
pr_pi1 = (primedec_t**)malloc(PRIME_MAXSIZE*sizeof(primedec_t*));
for(i=0; i<PRIME_MAXSIZE; i++) {
pr_pi1[i] = (primedec_t*)malloc(PRIME_MAXSIZE*sizeof(primedec_t));
}
q = me*ne - 1;
/* Fill-in pr_q, t, pr_pi1, t_pi1, gi, Nif, Kif and dif */
GKK_prepare(q, ne, pr_q, &t, pr_pi1, t_pi1, gi, Nif, Kif, dif);
#ifdef DEBUG_IPT
GKK_output_tables(me, ne, q, pr_q, t, gi, Nif, Kif, dif);
#endif
ht = t;
if( (q % 2) == 0 )
ht = t + 1;
ndiv = 1;
for(i=0; i<t; i++) {
fi[i] = pr_q[i].e;
pifi[i] = pr_q[i].pe;
ndiv = ndiv * (pr_q[i].e + 1);
}
/*
* First loop to compute number of leaders
*/
/* loop over divisors v of q */
nlead = 0;
for (div=0; div<ndiv-1; div++) {
GKK_L(t, pr_q, fi, Kif, dif, Li, diLi, &cl, &nl);
nlead += nl;
/* step to next divisor v of q (q/v = product of pifi) */
for(i=0; i<t; i++) {
if( fi[i] == 0 ) {
fi[i] = pr_q[i].e;
pifi[i] = pr_q[i].pe;
}
else {
fi[i] = fi[i] - 1;
pifi[i] = pifi[i] / pr_q[i].p;
break;
}
}
}
*nleaders = nlead;
/*fprintf(stderr, "Number of leader : %ld => ", (long)nlead); */
/*
* Second loop to store leaders
*/
(*leaders) = (int*) malloc (nlead*3*sizeof(int));
/* Restore fi and pifi */
for(i=0; i<t; i++) {
fi[i] = pr_q[i].e;
pifi[i] = pr_q[i].pe;
}
/* loop over divisors v of q */
nlead = 0;
for (div=0; div<ndiv-1; div++) {
GKK_L(t, pr_q, fi, Kif, dif, Li, diLi, &cl, &nl);
if( div == 0 ) {
GKK_precompute_terms(q, pr_q, t, gi, diLi, rp, Mg, SIZE_MG);
memcpy(vMg, Mg, SIZE_MG*sizeof(int));
}
else {
for(i=0; i<t; i++) {
if( pifi[i] != pr_q[i].pe ) {
for(j = rp[i]; j<rp[i]+diLi[i]; j++) {
vMg[j] = ((pr_q[i].pe / pifi[i])*Mg[j]) % q;
}
}
else {
memcpy(&(vMg[rp[i]]), &(Mg[rp[i]]), diLi[i]*sizeof(int));
}
}
}
for (i=0;i<ht; i++)
si[i] = 0;
while (1) {
/* convert current leader from additive to multiplicative domain */
s = 0;
for (i=0; i<t; i++) {
if( diLi[i] == 0 )
continue;
s = (s + vMg[rp[i] + si[i]]) % q;
}
/* assign it to owner */
(*leaders)[nlead ] = s;
(*leaders)[nlead+1] = cl;
//leaders[nleaders+2] = owner;
nlead += 3;
nl = nl - 1;
if( nl == 0 )
break;
/* step to next leader */
for(i=0; i<ht; i++) {
if( diLi[i] == 0 ) continue;
si[i] = si[i] + 1;
if( si[i] < diLi[i] ) {
if( (i == (ht-1)) && (ht != t) ) {
/* flip vMg for i=1 */
for(j=0; j<diLi[0]; j++) {
vMg[rp[0]+j] = (q-vMg[rp[0]+j]) % q;
}
}
break;
}
else
si[i] = 0;
}
}
/* step to next divisor v of q (q/v = product of pifi) */
for(i=0; i<t; i++) {
if( fi[i] == 0 ) {
fi[i] = pr_q[i].e;
pifi[i] = pr_q[i].pe;
}
else {
fi[i] = fi[i] - 1;
pifi[i] = pifi[i] / pr_q[i].p;
break;
}
}
}
#ifdef DEBUG_IPT
printf("nleaders : %d\n", *nleaders);
printf("leaders : ");
for(i=0; i<*nleaders; i++)
{
printf(" %d", (*leaders)[i*3]);
printf(" %d", (*leaders)[i*3+1]);
printf(" %d", (*leaders)[i*3+2]);
}
printf("\n");
#endif
for(i=0; i<PRIME_MAXSIZE; i++)
free(pr_pi1[i]);
free(pr_pi1);
free(ptr);
return PLASMA_SUCCESS;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void GKK_L ( int  t,
primedec_t *  pr_q,
int *  fi,
int *  Kif,
int *  dif,
int *  Li,
int *  diLi,
int *  cl,
int *  nl 
)

GKK_L computes Li, cycle length cl, and number of leaders nl.

See also:
GKK_prepare
Parameters:
[in]tNumber of primes in prime decomposition of q
[in]pr_qPrime Decomposition of q
[in]fiArray of dimension PWR_MAXSIZE
[in]KifArray of dimension PWR_MAXSIZE*PRIME_MAXSIZE Stores the multiplicative order of n modulo pi^ei
[in]difArray of dimension PWR_MAXSIZE*PRIME_MAXSIZE dif(i,j) = Nif(i,j) / Kif(i,j)
[out]LiArray of dimension PWR_MAXSIZE
[out]diLiArray of dimension PWR_MAXSIZE diLi(i) = di(i) * Li (i) (bound for the leaders sets)
[out]clCycle length
[out]nlNumber of leaders

Definition at line 462 of file gkkleader.c.

References gcd(), and PWR_MAXSIZE.

{
int i, lcl, lnl;
/* compute cycle length (cl) and Li */
lcl = 1;
for (i=0; i<t; i++) {
if( fi[i] == 0 ) {
Li[i] = 1;
}
else {
Li[i] = gcd( Kif[i*PWR_MAXSIZE+fi[i]-1], lcl);
lcl = (lcl * Kif[i*PWR_MAXSIZE+fi[i]-1]) / Li[i];
}
}
if( pr_q[0].p == 2 ) {
if( fi[0] == 0 ) {
Li[t] = 1;
}
else {
int tmp = Kif[t*PWR_MAXSIZE+fi[0]-1];
Li[t] = gcd( tmp, lcl);
lcl = (lcl * tmp) / Li[t];
}
}
*cl = lcl;
/* count number of cycles (nl) and compute di*Li */
lnl = 1;
for (i=0; i<t; i++) {
if( fi[i] != 0 ) {
diLi[i] = dif[i*PWR_MAXSIZE+fi[i]-1] * Li[i];
lnl = lnl * diLi[i];
}
else
diLi[i] = 0;
}
if( pr_q[0].p == 2 ) {
if( fi[0] != 0 ) {
diLi[t] = dif[t*PWR_MAXSIZE+fi[0]-1] * Li[t];
lnl = lnl * diLi[t];
}
else
diLi[t] = 0;
}
*nl = lnl;
#ifdef DEBUG_IPT
printf("nl : %d\n", lnl);
printf("cl : %d\n", lcl);
printf("Li :");
for(i=0; i<t+1; i++)
printf(" %d", Li[i]);
printf("\n");
printf("diLi :");
for(i=0; i<t+1; i++)
printf(" %d", diLi[i]);
printf("\n");
#endif
}

Here is the call graph for this function:

Here is the caller graph for this function:

int GKK_modpow ( int *  dt,
int  e,
int  m 
)

GKK_modpow computes x^e mod m, like the modpow function. Unlike modpow, x is given implicitly via its doubling table dt.

See also:
GKK_doublingtable
Parameters:
[in]dtArray of x^i mod m.
[in]e
[in]m
Returns:
Return values:
x^emod m

Definition at line 102 of file gkkleader.c.

{
int64_t y8, m8;
int z, i;
m8 = m;
y8 = 1;
z = e;
i = 0;
/* dt size is not tested because tested before in GKK_doublingtable */
while( z > 0 ) {
if( (z%2) == 1 )
y8 = (y8*dt[i]) % m8;
z = z / 2;
i++;
}
return (int)y8;
}

Here is the caller graph for this function:

int GKK_multorder ( int  n,
int  p,
int  e,
int  pe,
primedec_t *  pr_p1,
int  t_p1 
)

GKK_multorder finds the multiplicative order of n modulo p^e = pe, given the prime decomposition of p-1.

See also:
GKK_doublingtable
GKK_modpow
Parameters:
[in]n
[in]p
[in]e
[in]pe
[in,out]pr_p1On entry, prime decomposition of p - 1 On exit, prime decomposition of p - 1 and p^(e-1) appends at the end
[in]t_p1Number of primes in prime decomposition of p-1
Returns:
Return values:
themultiplicative order of n modulo p^e = pe

Definition at line 232 of file gkkleader.c.

References primedec::e, GKK_doublingtable(), GKK_modpow(), primedec::p, primedec::pe, PRIME_MAXSIZE, and PWR_MAXSIZE.

{
int fi[PRIME_MAXSIZE], dt[PWR_MAXSIZE];
int t, i, K1, c;
int K;
/* append p^(e-1) to the prime decomposition of p-1 */
t = t_p1;
if( e > 1 ) {
pr_p1[t].p = p;
pr_p1[t].e = e-1;
pr_p1[t].pe = pe / t;
t = t_p1 + 1;
}
GKK_doublingtable(n, pe, pe-1, dt);
for (i=0; i<PRIME_MAXSIZE; i++)
fi[i] = pr_p1[i].e;
K = pe - pe/p;
for(i=0; i<t; i++) {
while( fi[i] > 0 ) {
K1 = K / pr_p1[i].p;
c = GKK_modpow(dt, K1, pe);
if( c == 1 ) {
K = K1;
fi[i] = fi[i] - 1;
}
else {
break;
}
}
}
return K;
}

Here is the call graph for this function:

Here is the caller graph for this function:

void GKK_output_tables ( int  m,
int  n,
int  q,
primedec_t *  pr_q,
int  t,
int *  gi,
int *  Nif,
int *  Kif,
int *  dif 
)

GKK_output_tables prints information from the heart of the GKK algorithm in a human friendly format.

Parameters:
[in]m
[in]n
[in]q
[in]pr_qPrime Decomposition of q
[in]tNumber of primes in prime decomposition of q
[in]giArray of size PRIME_MAXSIZE Stores prime root of each pi^ei
[in]NifArray of dimension PWR_MAXSIZE*PRIME_MAXSIZE Stores the Phi(pi^fi) with fi in (1 .. ei)
[in]KifArray of dimension PWR_MAXSIZE*PRIME_MAXSIZE Stores the multiplicative order of n modulo pi^ei
[in]difArray of dimension PWR_MAXSIZE*PRIME_MAXSIZE dif(i,j) = Nif(i,j) / Kif(i,j)

Definition at line 707 of file gkkleader.c.

References primedec::e, lapack_testing::f, and PWR_MAXSIZE.

{
int i, f;
fprintf(stdout, "Information from the GKK algorithm\n");
fprintf(stdout, "==================================\n");
fprintf(stdout, "m = %4d\n", m);
fprintf(stdout, "n = %4d\n", n);
fprintf(stdout, "q = %4d\n", q);
for (i=0; i<t; i++) {
fprintf(stdout, "==================================\n");
fprintf(stdout, " i = %4d\n", i+1 );
fprintf(stdout, " p = %4d\n", pr_q[i].p);
fprintf(stdout, " e = %4d\n", pr_q[i].e);
fprintf(stdout, " p^e = %4d\n", pr_q[i].pe);
if( pr_q[i].p != 2 )
fprintf(stdout, " g = %4d\n", gi[i]);
else
fprintf(stdout, "mod(n,4) = %4d\n", n%4);
fprintf(stdout, "\n");
fprintf(stdout, " f | ");
for (f=0; f<pr_q[i].e; f++)
fprintf(stdout, "%4d ", f+1);
fprintf(stdout, "\n");
fprintf(stdout, "----------------------------------\n");
fprintf(stdout, "Ni(f) | ");
for (f=0; f<pr_q[i].e; f++)
fprintf(stdout, "%4d ", Nif[i*PWR_MAXSIZE+f]);
fprintf(stdout, "\n");
fprintf(stdout, "Ki(f) | ");
for (f=0; f<pr_q[i].e; f++)
fprintf(stdout, "%4d ", Kif[i*PWR_MAXSIZE+f]);
fprintf(stdout, "\n");
fprintf(stdout, "di(f) | ");
for (f=0; f<pr_q[i].e; f++)
fprintf(stdout, "%4d ", dif[i*PWR_MAXSIZE+f]);
fprintf(stdout, "\n");
fprintf(stdout, "\n");
}
}

Here is the caller graph for this function:

void GKK_precompute_terms ( int  q,
primedec_t *  pr_q,
int  t,
int *  gi,
int *  diLi,
int *  rp,
int *  Mg,
int  nMg 
)

GKK_precompute_terms Precompute M*g^i mod q.

See also:
GKK_prepare
Parameters:
[in]q
[in]pr_qPrime Decomposition of q
[in]tNumber of primes in prime decomposition of q
[in]giArray of dimension PWR_MAXSIZE
[in]diLiArray of dimension PWR_MAXSIZE
[out]rpArray of dimension PWR_MAXSIZE rp(i) = sum(diLi(i), i=0..i-1), rp(0) = 0;
[out]MgArray of dimension nMg
[in]nMg

Definition at line 563 of file gkkleader.c.

References primedec::pe, and plasma_error().

{
int i, j, nMg_needed, ht;
ht = t;
if( pr_q[0].p == 2 )
ht = t + 1;
nMg_needed = diLi[0];
for(i=1; i<ht; i++)
nMg_needed += diLi[i];
if( nMg_needed > nMg ) {
plasma_error(__func__, "the size of Mg is not large enough");
return;
}
rp[0] = 0;
for(i=0; i<t; i++) {
if( pr_q[i].p != 2 ) {
Mg[rp[i]] = q / pr_q[i].pe;
for(j=1; j<diLi[i]; j++)
Mg[rp[i]+j] = (Mg[rp[i]+j-1]*gi[i]) % q;
}
rp[i+1] = rp[i] + diLi[i];
}
if( pr_q[0].p == 2 ) {
Mg[rp[0]] = q / pr_q[0].pe;
for(j=1; j<diLi[0]; j++)
Mg[rp[0]+j] = (Mg[rp[0]+j-1]*5) % q;
}
#ifdef DEBUG_IPT
printf("t : %d\n", t);
printf("ht : %d\n", ht);
printf("rp :");
for(i=0; i<ht; i++)
printf(" %d", rp[i]);
printf("\n");
printf("Mg :");
for(i=0; i<nMg; i++)
printf(" %d", Mg[i]);
printf("\n");
#endif
}

Here is the call graph for this function:

Here is the caller graph for this function:

void GKK_prepare ( int  q,
int  n,
primedec_t *  pr_q,
int *  t,
primedec_t **  pr_pi1,
int *  t_pi1,
int *  gi,
int *  Nif,
int *  Kif,
int *  dif 
)

GKK_prepare Prepare for the generation of leaders.

See also:
GKK_primroot
GKK_multorder
factor
gcd
Parameters:
[in]q= mn - 1
[in]n
[out]pr_qPrime Decomposition of q
[out]tNumber of primes in prime decomposition of q
[out]pr_pi1Array of prime decompositions of each pi in prime decomposition of q
[out]t_pi1Array of number of primes in each prime decomposition of pi found in prime decomposition of q
[out]giArray of size PRIME_MAXSIZE Stores prime root of each pi^ei
[out]NifArray of dimension PWR_MAXSIZE*PRIME_MAXSIZE Stores the Phi(pi^fi) with fi in (1 .. ei) Phi(k) is the number of positive integers less than or equal to k which are coprime to k.
[out]KifArray of dimension PWR_MAXSIZE*PRIME_MAXSIZE Stores the multiplicative order of n modulo pi^ei
[out]difArray of dimension PWR_MAXSIZE*PRIME_MAXSIZE dif(i,j) = Nif(i,j) / Kif(i,j)

Definition at line 318 of file gkkleader.c.

References primedec::e, lapack_testing::f, factor(), gcd(), GKK_multorder(), GKK_primroot(), primedec::p, and PWR_MAXSIZE.

{
int i, f, p1f, tmp;
/* factor q */
factor(q, pr_q, t);
/* factor pi-1 for each pi in q */
for (i=0; i<*t; i++)
factor(pr_q[i].p-1, pr_pi1[i], &(t_pi1[i]));
/*
* Compute the number of positive integers coprime
* to f, Ni(f) for f = 1, 2, ..., ei
*/
for (i=0; i<*t; i++) {
tmp = PWR_MAXSIZE*i;
Nif[tmp] = pr_q[i].p - 1;
for(f=1; f < pr_q[i].e; f++)
Nif[tmp+f] = pr_q[i].p * Nif[tmp+f-1];
}
if( pr_q[0].p == 2 ) {
for(f=1; f < pr_q[0].e; f++)
Nif[f] = Nif[f] / 2;
Nif[PWR_MAXSIZE*(*t)] = 1;
for(f=1; f < pr_q[0].e; f++)
Nif[PWR_MAXSIZE*(*t)+f] = 2;
}
/*
* Case A: odd pi
*/
for(i=0; i <*t; i++) {
if( pr_q[i].p != 2 ) {
tmp = PWR_MAXSIZE*i;
/* compute Ki(ei) */
Kif[tmp+pr_q[i].e-1] = GKK_multorder(n, pr_q[i].p, pr_q[i].e, pr_q[i].pe, pr_pi1[i], t_pi1[i]);
/* compute Ki(f) for f = 1, 2, ..., ei-1 */
for(f = pr_q[i].e-2; f>-1; f--) {
if( Kif[tmp+f+1] >= pr_q[i].p )
Kif[tmp+f] = Kif[tmp+f+1] / pr_q[i].p;
else
Kif[tmp+f] = Kif[tmp+f+1];
}
/* compute di(f) for f = 1, 2, ..., ei */
for(f=0; f<pr_q[i].e; f++)
dif[tmp+f] = Nif[tmp+f] / Kif[tmp+f];
/* compute primitive root gi for pi^ei */
if( dif[tmp+pr_q[i].e-1] == 1 )
/* special case where n is already a primitive root*/
gi[i] = n % pr_q[i].pe;
else
/* general case requiring a search for a primitive root */
gi[i] = GKK_primroot(pr_q[i].p, pr_q[i].e, pr_pi1[i], t_pi1[i]);
}
}
/*
* Case B: even pi (only the first one)
*/
if( pr_q[0].p == 2 ) {
/* reset primitive root */
gi[0] = 0;
/* compute d1(f) for f = 1,2,...,e1 */
p1f = 2;
for (f=0; f<pr_q[0].e; f++) {
if( (n%4) == 1)
dif[f] = gcd( ( (n%p1f)-1)/4, Nif[f] );
else
dif[f] = gcd( (p1f-(n%p1f)-1)/4, Nif[f] );
p1f = p1f*2;
}
/* compute K1(f) for f = 1,2,...,e1 */
for (f=0; f<pr_q[0].e; f++)
Kif[f] = Nif[f] / dif[f];
/* compute K0(f) for f = 1,2,...,e1 */
tmp = PWR_MAXSIZE * (*t);
Kif[tmp] = 1;
for (f=1; f<pr_q[0].e; f++) {
Kif[tmp+f] = 2;
if( (n%4) == 1 )
Kif[tmp+f] = 1;
}
/* compute d0(f) for f = 1,2,...,e1 */
for (f=0; f<pr_q[0].e; f++)
dif[tmp+f] = Nif[tmp+f] / Kif[tmp+f];
}
}

Here is the call graph for this function:

Here is the caller graph for this function:

int GKK_primroot ( int  p,
int  e,
primedec_t *  pr_p1,
int  t_p1 
)

Definition at line 150 of file gkkleader.c.

References GCAND_SIZE, GKK_doublingtable(), GKK_modpow(), plasma_error(), and PWR_MAXSIZE.

{
static int gcand[GCAND_SIZE] = {
2, 3, 5, 6, 7, 10, 11, 12, 13, 14,
15, 17, 18, 19, 20, 21, 22, 23, 24, 26,
28, 29, 30, 31, 33, 34, 35, 37, 38, 39,
41, 43, 44, 46, 47, 51, 52, 53, 55, 57,
58, 59, 60, 62, 69, 73 };
int dt[PWR_MAXSIZE];
int i, j, p2, c, g;
int ok;
p2 = p;
if( e >= 2 )
p2 = p*p;
/* try candidates */
for(i=0; i < GCAND_SIZE; i++) {
g = gcand[i];
GKK_doublingtable(g, p2, p-1, dt);
ok = 1;
for(j=0; j<t_p1; j++) {
c = GKK_modpow(dt, (p-1) / pr_p1[j].p, p2);
if( (c % p) == 1 ) {
ok = 0;
break;
}
}
if( ok == 1 )
break;
}
if( ok != 1 ) {
plasma_error(__func__, "failed to find a primitive root");
return -1;
}
/* check that g is a primitive root also for p^k for all k */
if( e >= 2 ) {
c = GKK_modpow(dt, p-1, p2);
if( c == 1 ) {
/* g is not a primitive root for p^k, but g+p is */
g = g + p;
}
}
return g;
}

Here is the call graph for this function:

Here is the caller graph for this function: