1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91
|
#define Nbits 40
unsigned long long L1, L2;
unsigned short W1[Nbits], W2[Nbits],
W1_sqr[2*Nbits], W2_sqr[2*Nbits],
W_sm[2*(Nbits+1)]; // en principe 2*Nbits+1 suffisent mais un nombre paire est + facile a programmer dans 'pendule'
void Bit40( unsigned long long L, unsigned short * W)
{
short i;
for (i=0; i < Nbits; i++)
{
W[i] = L & 1;
L = L>> 1;
}
}
void LSQR( unsigned short * W,unsigned short * W2)
{
short i;
memset(W2,0,2*Nbits*sizeof(short));
for (i=0; i< Nbits; i++) if ( W[i] )
{
short retenue =0;
short k;
for (k=0; k < Nbits; k++)
{
short u= W2[i+k] + W[k] + retenue;
if ( u <= 1)
{
W2[i+k] = u;
retenue = 0;
}
else
{
W2[i+k] = u - 2;
retenue = 1;
}
}
}
}
void Add2( unsigned short * W1_sqr, unsigned short * W2_sqr,unsigned short * W_sm)
{
short k, retenue = 0;
memset(W_sm,0,2*(Nbits+1)*sizeof(short));
for (k=0; k < 2*Nbits; k++)
{
short m = retenue+W1_sqr[k] + W2_sqr[k] ;
if ( m <=1)
{
W_sm[k] = m;
retenue = 0;
}
else
{
W_sm[k] = m-2;
retenue = 1;
}
}
W_sm[2*Nbits]=retenue; //W_sm[2*Nbits+1] est toujours nul
}
void potence (unsigned short *W_sm, unsigned long long * square, unsigned long long * reste)
{
*square=0;
* reste =0;
short i = Nbits;
while (1)
{
long long L;
short j=2*i+1;
short k=2*i;
* reste = (* reste) * 4 + (W_sm[j] << 1) + W_sm[k];
L = 4* (*square) + 1;
if ( L <= (* reste) )
{
(* reste) -= L;
(*square) = 2*(*square) +1;
}
else
(* square) = 2 * ( * square);
i--;
if ( i < 0) break;
}
}
void Norme(unsigned long long L1, unsigned long long L2, unsigned long long * L, unsigned long long * R)
{
Bit40(L1,W1);
Bit40(L2,W2);
LSQR(W1,W1_sqr);
LSQR(W2,W2_sqr);
Add2(W1_sqr,W2_sqr,W_sm);
potence (W_sm, L,R);
} |
Partager