libtommath/bn_mp_compute_factored_factorial.c at master · czurnieden/libtommath · GitHub

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
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
#include <tommath.h>
#ifdef BN_MP_COMPUTE_FACTORED_FACTORIAL_C

static long local_highbit(unsigned  long n)
{
   long r=0;
   while (n >>= 1) {
      r++;
   }
   return r;
}

unsigned long BN_MP_FACTORED_FACTORIAL_CUTOFF = 0x10000000;
/* Input is to be sorted as: prime,exponent,prime,exponent,... */
int mp_compute_factored_factorial(unsigned long *f, unsigned long f_length,
                                  mp_int *c,   unsigned long stop)
{

   unsigned long length,start,i;
   unsigned long shift = 0;
   long bit,k=0,hb;
   mp_int temp;
   int e;

   if (stop == 0) {
      length = f_length;
   } else {
      length = stop;
   }
   if (f[0] == 2 && f[1] > 0) {
      shift = f[1];
      if (f_length == 2) {
         // TODO: goto final_shift;
         if ((e = mp_set_int(c,1)) != MP_OKAY) {
            return e;
         }
         if ((e = mp_mul_2d(c,(int)f[1],c)) != MP_OKAY) {
            return e;
         }
         return MP_OKAY;
      }
   }

   start = 0;
   if (shift) {
      start+=2;
      k=2;
   }
   bit = 0;
   for (; k<(long)f_length; k+=2) {
      hb=local_highbit(f[k+1]);
      if (bit < hb)bit=hb;
   }
   /* just for safety reasons, y'know */
   if (bit >(long) DIGIT_BIT) {
      return MP_VAL;
   }
   if ((e = mp_set_int(c,1)) != MP_OKAY) {
      return e;
   }
   if (f_length > BN_MP_FACTORED_FACTORIAL_CUTOFF) {
      if ((e = mp_init(&temp)) != MP_OKAY) {
         return e;
      }
      for (; bit>=0; bit--) {
         if ((e =  mp_sqr(c, c)) != MP_OKAY) {
            return e;
         }
         if ((e =  mp_set_int(&temp,1)) != MP_OKAY) {
            return e;
         }
         for (i=start; i<f_length; i+=2) {
            if ((f[i+1] & (1LU << (unsigned long)bit)) != 0) {
               if ((e =  mp_mul_d(&temp,(mp_digit)f[i], &temp)) != MP_OKAY) {
                  return e;
               }
            }
         }
         if ((e =  mp_mul(&temp,c,c)) != MP_OKAY) {
            return e;
         }
      }
   } else {
      for (; bit>=0; bit--) {
         if ((e = mp_sqr(c, c)) != MP_OKAY) {
            return e;
         }
         for (i=start; i<f_length; i+=2) {
            if ((f[i+1] & (1LU << (unsigned long)bit)) != 0) {
               if ((e = mp_mul_d(c,(mp_digit)f[i], c)) != MP_OKAY) {
                  return e;
               }
            }
         }
      }
   }
// TODO: final_shift:
   if (shift) {
      /* DIGIT_BIT might be 60 which makes this extra check necessary */
      if (shift < INT_MAX) {
         /* The choice of types is a bit questionable. Sometimes. */
         if ((e = mp_mul_2d(c, (int) shift, c)) != MP_OKAY) {
            return e;
         }
      } else {
         int multiplicator=0;
         while (shift > INT_MAX) {
            shift >>= 1;
            multiplicator++;
         }
         if ((e = mp_mul_2d(c, (int) shift, c)) != MP_OKAY) {
            return e;
         }
         if ((e = mp_mul_2d(c, 1<<multiplicator, c)) != MP_OKAY) {
            return e;
         }
      }
   }
   return MP_OKAY;
}
#endif