source: git/factory/cf_binom.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: 3.3 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id$ */
3
4#include "config.h"
5
6#include "cf_assert.h"
7
8#include "cf_defs.h"
9#include "canonicalform.h"
10#include "cf_binom.h"
11
12#define MAXPT 40
13
14#define INITPT 10
15
16CFArray * ptZ = 0;
17CFArray * ptF = 0;
18
19int charac = 0;
20int gfdeg = 0;
21int ptZmax = INITPT;
22int ptFmax = 0;
23
24void
25resetFPT()
26{
27    ptFmax = 0;
28}
29
30void
31initPT ( )
32{
33    static bool initialized = false;
34
35    if ( ! initialized ) {
36        initialized = true;
37        ptZ = new CFArray[MAXPT+1];
38        ptF = new CFArray[MAXPT+1];
39        int i, j;
40        ptZ[0] = CFArray(1); ptZ[0][0] = 1;
41        ptF[0] = CFArray(1);
42        for ( i = 1; i <= INITPT; i++ ) {
43            ptF[i] = CFArray(i+1);
44            ptZ[i] = CFArray(i+1);
45            (ptZ[i])[0] = 1;
46            for ( j = 1; j < i; j++ )
47                (ptZ[i])[j] = (ptZ[i-1])[j-1] + (ptZ[i-1])[j];
48            (ptZ[i])[i] = 1;
49        }
50        for ( i = INITPT+1; i <= MAXPT; i++ ) {
51            ptF[i] = CFArray(i+1);
52            ptZ[i] = CFArray(i+1);
53        }
54        ptZmax = INITPT;
55        ptFmax = 0;
56    }
57}
58
59CanonicalForm
60binomialpower ( const Variable & x, const CanonicalForm & a, int n )
61{
62    if ( n == 0 )
63        return 1;
64    else if ( n == 1 )
65        return x + a;
66    else if ( getCharacteristic() == 0 ) {
67        if ( n <= MAXPT ) {
68            if ( n > ptZmax ) {
69                int i, j;
70                for ( i = ptZmax+1; i <= n; i++ ) {
71                    (ptZ[i])[0] = 1;
72                    for ( j = 1; j < i; j++ )
73                        (ptZ[i])[j] = (ptZ[i-1])[j-1] + (ptZ[i-1])[j];
74                    (ptZ[i])[i] = 1;
75                }
76                ptZmax = n;
77            }
78            CanonicalForm result = 0, apower = 1;
79            int k;
80            for ( k = n; k >= 0; k-- ) {
81                result += power( x, k ) * apower * (ptZ[n])[k];
82                if ( k != 0 )
83                    apower *= a;
84            }
85            return result;
86        }
87        else {
88            CanonicalForm result = binomialpower( x, a, MAXPT );
89            CanonicalForm xa = x + a;
90            int i;
91            for ( i = MAXPT; i < n; i++ )
92                result *= xa;
93            return result;
94        }
95    }
96    else {
97        if ( getCharacteristic() != charac || gfdeg != getGFDegree() ) {
98            ptFmax = 0;
99            charac = getCharacteristic();
100            gfdeg = getGFDegree();
101            (ptF[0])[0] = 1;
102        }
103        if ( n <= MAXPT ) {
104            if ( n > ptFmax ) {
105                int i, j;
106                for ( i = ptFmax+1; i <= n; i++ ) {
107                    (ptF[i])[0] = 1;
108                    for ( j = 1; j < i; j++ )
109                        (ptF[i])[j] = (ptF[i-1])[j-1] + (ptF[i-1])[j];
110                    (ptF[i])[i] = 1;
111                }
112                ptFmax = n;
113            }
114            CanonicalForm result = 0, apower = 1;
115            int k;
116            for ( k = n; k >= 0; k-- ) {
117                result += power( x, k ) * apower * (ptF[n])[k];
118                if ( k != 0 )
119                    apower *= a;
120            }
121            return result;
122        }
123        else {
124            CanonicalForm result = binomialpower( x, a, MAXPT );
125            CanonicalForm xa = x + a;
126            int i;
127            for ( i = MAXPT; i < n; i++ )
128                result *= xa;
129            return result;
130        }
131    }
132}
Note: See TracBrowser for help on using the repository browser.