source: git/factory/ffops.cc @ fd80670

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