source: git/factory/ffops.cc @ 48e0bcb

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