source: git/factory/ffops.cc @ cd86ac

spielwiese
Last change on this file since cd86ac was a43bd1, checked in by Olaf Bachmann <obachman@…>, 24 years ago
2000-08-04 Olaf Bachmann <obachman@mathematik.uni-kl.de> * ffops.cc: initialize ff_prime to 0, so that first call to ff_setprime always initializes ff_invtab. git-svn-id: file:///usr/local/Singular/svn/trunk@4507 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 2.1 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id: ffops.cc,v 1.9 2000-08-04 11:10:15 obachman Exp $ */
3
4#include <config.h>
5
6#include <string.h>
7
8#include "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
18#ifdef __MWERKS__
19#include <stuff_64.h>
20int ff_mul( const int a, const int b )
21{
22  Int_64 c;
23  unsigned int ua, ub;
24  int res;
25 
26  if ((!a) || (!b)) return 0;
27  if (a>0)
28    ua = a;
29  else
30    ua = -a;
31  if (b>0)
32    ub = b;
33  else
34    ub = -b;
35  mul_64(&c, ua, ub);
36  res = rem_64(&c, ff_prime);
37  if ((a>0) != (b>0))
38    res = ff_prime-res;
39  return res;
40}
41#endif
42
43void ff_setprime ( const int p )
44{
45    if ( p != ff_prime ) {
46        ff_prime = p;
47        ff_halfprime = ff_prime / 2;
48        if ( ! ff_big )
49            memset( ff_invtab, 0, ff_prime*sizeof(short) );
50    }
51}
52
53int ff_newinv ( const int a )
54{
55    int p, q, r1, r2, y1, y2;
56    if (a < 2)
57          return (ff_invtab[a] = a);
58    r1 = p = ff_prime;
59    q = r1 / a;
60    y1 = -q;
61    r1 -= a * q;
62    if (r1 == 1)
63    {
64      y1 += p;
65      ff_invtab[y1] = a;
66      return (ff_invtab[a] = y1);
67    }
68    r2 = a;
69    y2 = 1;
70    for (;;)
71    {
72      q = r2 / r1;
73      y2 -= y1 * q;
74      r2 -= r1 * q;
75      if (r2 == 1)
76      {
77        if (y2 < 0)
78          y2 += p;
79        ff_invtab[y2] = a;
80        return (ff_invtab[a] = y2);
81      }
82      q = r1 / r2;
83      y1 -= y2 * q;
84      r1 -= r2 * q;
85      if (r1 == 1)
86      {
87        if (y1 < 0)
88          y1 += p;
89        ff_invtab[y1] = a;
90        return (ff_invtab[a] = y1);
91      }
92    }
93}
94
95int ff_biginv ( const int a )
96{
97    int p, q, r1, r2, y1, y2;
98    if (a < 2)
99      return a;
100    r1 = p = ff_prime;
101    q = r1 / a;
102    y1 = -q;
103    r1 -= a * q;
104    if (r1 == 1)
105      return p + y1;
106    r2 = a;
107    y2 = 1;
108    for (;;)
109    {
110      q = r2 / r1;
111      y2 -= y1 * q;
112      r2 -= r1 * q;
113      if (r2 == 1)
114      {
115        if (y2 > 0)
116          return y2;
117        else
118          return p + y2;
119      }
120      q = r1 / r2;
121      y1 -= y2 * q;
122      r1 -= r2 * q;
123      if (r1 == 1)
124      {
125        if (y1 > 0)
126          return y1;
127        else
128          return p + y1;
129      }
130    }
131}
Note: See TracBrowser for help on using the repository browser.