source: git/factory/ffops.cc @ e4fe2b

jengelh-datetimespielwiese
Last change on this file since e4fe2b was e4fe2b, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
FIX: Fixed huge BUG in cf_gmp.h CHG: starting to cleanup factory
  • Property mode set to 100644
File size: 1.8 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id$ */
3
4#include "config.h"
5
6#include <string.h>
7
8#include "cf_assert.h"
9
10#include "cf_defs.h"
11#include "ffops.h"
12
13int ff_prime = 0;
14int ff_halfprime = 0;
15bool ff_big = false;
16short * ff_invtab = new short [32767];
17
18void ff_setprime ( const int p )
19{
20    if ( p != ff_prime ) {
21        ff_prime = p;
22        ff_halfprime = ff_prime / 2;
23        if ( ! ff_big )
24            memset( ff_invtab, 0, ff_prime*sizeof(short) );
25    }
26}
27
28int ff_newinv ( const int a )
29{
30    int p, q, r1, r2, y1, y2;
31    if (a < 2)
32          return (ff_invtab[a] = a);
33    r1 = p = ff_prime;
34    q = r1 / a;
35    y1 = -q;
36    r1 -= a * q;
37    if (r1 == 1)
38    {
39      y1 += p;
40      ff_invtab[y1] = a;
41      return (ff_invtab[a] = y1);
42    }
43    r2 = a;
44    y2 = 1;
45    for (;;)
46    {
47      q = r2 / r1;
48      y2 -= y1 * q;
49      r2 -= r1 * q;
50      if (r2 == 1)
51      {
52        if (y2 < 0)
53          y2 += p;
54        ff_invtab[y2] = a;
55        return (ff_invtab[a] = y2);
56      }
57      q = r1 / r2;
58      y1 -= y2 * q;
59      r1 -= r2 * q;
60      if (r1 == 1)
61      {
62        if (y1 < 0)
63          y1 += p;
64        ff_invtab[y1] = a;
65        return (ff_invtab[a] = y1);
66      }
67    }
68}
69
70int ff_biginv ( const int a )
71{
72    int p, q, r1, r2, y1, y2;
73    if (a < 2)
74      return a;
75    r1 = p = ff_prime;
76    q = r1 / a;
77    y1 = -q;
78    r1 -= a * q;
79    if (r1 == 1)
80      return p + y1;
81    r2 = a;
82    y2 = 1;
83    for (;;)
84    {
85      q = r2 / r1;
86      y2 -= y1 * q;
87      r2 -= r1 * q;
88      if (r2 == 1)
89      {
90        if (y2 > 0)
91          return y2;
92        else
93          return p + y2;
94      }
95      q = r1 / r2;
96      y1 -= y2 * q;
97      r1 -= r2 * q;
98      if (r1 == 1)
99      {
100        if (y1 > 0)
101          return y1;
102        else
103          return p + y1;
104      }
105    }
106}
Note: See TracBrowser for help on using the repository browser.