File: mp_poly.C

package info (click to toggle)
gap-float 0.9.1%2Bds-6
  • links: PTS, VCS
  • area: main
  • in suites: bookworm, bullseye
  • size: 612 kB
  • sloc: ansic: 2,537; cpp: 1,998; xml: 184; makefile: 103; sh: 1
file content (250 lines) | stat: -rw-r--r-- 8,100 bytes parent folder | download | duplicates (4)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
/****************************************************************************
 *
 * mp_poly.C                                                Laurent Bartholdi
 *
 *   @(#)$id: fr_dll.c,v 1.18 2010/10/26 05:19:40 gap exp $
 *
 * Copyright (c) 2011, Laurent Bartholdi
 *
 ****************************************************************************
 *
 * driver for cpoly.C, using mpc_t complex numbers
 *
 ****************************************************************************/
#include <mpc.h>
#include <limits.h>
#include <math.h>

#define MPFR_REALS

#define cpoly cpoly_MPC

static mp_prec_t default_prec = 128; // ugly, non-reentrant

#ifdef MPFR_REALS
struct xreal {
  mpfr_t z;

  static const mp_rnd_t default_rnd;
  static const mp_prec_t default_prec = 32;

  // constructor
#ifdef DEBUG_MALLOC
  xreal() { printf("alloc() %p...",z); mpfr_init2(z,default_prec); printf("done\n"); }
  xreal(const double a){ printf("alloc(double %lf) %p...",a,z); mpfr_init2(z,default_prec); mpfr_set_d(z,a,default_rnd); printf("done\n"); }
  xreal(int m, int e){ printf("alloc(int %d,int %d) %p...",m,e,z); mpfr_init2(z,default_prec); mpfr_set_si_2exp(z,m,e,default_rnd); printf("done\n"); }
  xreal(const xreal &newz) { printf("alloc(xreal %p) %p...",newz.z,z); mpfr_init2(z,default_prec); mpfr_set_fr (z,newz.z,default_rnd); printf("done\n"); }
  ~xreal() { printf("free() %p...",z); mpfr_clear(z); printf("done\n"); }
#else
  xreal() { mpfr_init2(z,default_prec); }
  xreal(const double a){ mpfr_init2(z,default_prec); mpfr_set_d(z,a,default_rnd); }
  xreal(int m, int e){ mpfr_init2(z,default_prec); mpfr_set_si_2exp(z,m,e,default_rnd); }
  xreal(const xreal &newz) { mpfr_init2(z,default_prec); mpfr_set_fr (z,newz.z,default_rnd); }
  ~xreal() { mpfr_clear(z); }
#endif

  // operations
  xreal operator - () const { xreal newz; mpfr_neg(newz.z,z,default_rnd); return newz; };

  void operator += (const xreal &a){ mpfr_add(z,z,a.z,default_rnd); };
  void operator -= (const xreal &a){ mpfr_sub(z,z,a.z,default_rnd); };
  void operator *= (const xreal &a){ mpfr_mul(z,z,a.z,default_rnd); };
  void operator /= (const xreal &a){ mpfr_div(z,z,a.z,default_rnd); };
  
  void operator = (const mpfr_t &newz){ mpfr_set_prec(z,mpfr_get_prec(newz)); mpfr_set(z,newz,default_rnd) ;
  };
  void operator = (const xreal &newz){ mpfr_set_prec(z,mpfr_get_prec(newz.z)); mpfr_set(z,newz.z,default_rnd); };
};

bool const operator == (const xreal &a, const xreal &b) {
  return(mpfr_cmp(a.z,b.z) == 0);
}
bool const operator < (const xreal &a, const xreal &b) {
  return(mpfr_cmp(a.z,b.z) < 0);
}
bool const operator > (const xreal &a, const xreal &b) {
  return(mpfr_cmp(a.z,b.z) > 0);
}
bool const operator >= (const xreal &a, const xreal &b) {
  return(mpfr_cmp(a.z,b.z) >= 0);
}
bool const operator <= (const xreal &a, const xreal &b) {
  return(mpfr_cmp(a.z,b.z) <= 0);
}
bool const operator != (const xreal &a, const xreal &b) {
  return(mpfr_cmp(a.z,b.z) != 0);
}
xreal const operator + (const xreal &a, const xreal &b) {
  xreal newz; mpfr_add(newz.z,a.z,b.z,xreal::default_rnd); return newz;
}
xreal const operator - (const xreal &a, const xreal &b) {
  xreal newz; mpfr_sub(newz.z,a.z,b.z,xreal::default_rnd); return newz;
}
xreal const operator * (const xreal &a, const xreal &b) {
  xreal newz; mpfr_mul(newz.z,a.z,b.z,xreal::default_rnd); return newz;
}
xreal const operator / (const xreal &a, const xreal &b) {
  xreal newz; mpfr_div(newz.z,a.z,b.z,xreal::default_rnd); return newz;
}
#else
typedef double xreal;
#endif

struct xcomplex {
  mpc_t z;

  static const mp_rnd_t default_rnd;

  // constructor
  xcomplex(){ mpc_init2(z,default_prec); }
  xcomplex(const int i){ mpc_init2(z,default_prec); mpc_set_si(z,i,default_rnd); }
  xcomplex(const xcomplex &x){ mpc_init2(z,default_prec); mpc_set(z,x.z,default_rnd); }
#ifdef MPFR_REALS
  xcomplex(const double a){ mpc_init2(z,default_prec); mpc_set_d(z,a,default_rnd); }
  xcomplex(const xreal r){ mpc_init2(z,default_prec); mpc_set_fr(z,r.z,default_rnd); };
  xcomplex(const xreal a,const xreal b){ mpc_init2(z,default_prec); mpc_set_fr_fr(z,a.z,b.z,default_rnd); };
#else
  xcomplex(const xreal r){ mpc_init2(z,default_prec); mpc_set_d(z,r,default_rnd); }
  xcomplex(const xreal a,const xreal b){ mpc_init2(z,default_prec); mpc_set_d_d(z,a,b,default_rnd); }
#endif
  ~xcomplex() { mpc_clear(z); }

  // operations
  xcomplex operator - () const { xcomplex newz; mpc_neg(newz.z,z,default_rnd); return newz; };

  void operator += (const xcomplex &a){ mpc_add(z,z,a.z,default_rnd); };
  void operator -= (const xcomplex &a){ mpc_sub(z,z,a.z,default_rnd); };
  void operator *= (const xcomplex &a){ mpc_mul(z,z,a.z,default_rnd); };
  void operator /= (const xcomplex &a){ mpc_div(z,z,a.z,default_rnd); };
  
  void operator = (const mpc_t &newz){ mpc_set_prec(z,mpc_get_prec(newz)); mpc_set(z,newz,default_rnd) ;
  };
  void operator = (const xcomplex &newz){ mpc_set_prec(z,mpc_get_prec(newz.z)); mpc_set(z,newz.z,default_rnd); };
};

const mp_rnd_t xcomplex::default_rnd = mpfr_get_default_rounding_mode();
const long int xMAX_EXP = mpfr_get_emax();
const long int xMIN_EXP = mpfr_get_emin();
#ifdef MPFR_REALS
const mp_rnd_t xreal::default_rnd = mpfr_get_default_rounding_mode();
const xreal xINFIN(1,xMAX_EXP);
#else
const xreal xINFIN = pow(2,xMAX_EXP);
#endif

bool const operator == (const xcomplex &a, const xcomplex &b) {
  return(mpc_cmp(a.z,b.z) == 0);
}
bool const operator != (const xcomplex &a, const xcomplex &b) {
  return(mpc_cmp(a.z,b.z) != 0);
}
xcomplex const operator + (const xcomplex &a, const xcomplex &b) {
  xcomplex newz; mpc_add(newz.z,a.z,b.z,xcomplex::default_rnd); return newz;
}
xcomplex const operator - (const xcomplex &a, const xcomplex &b) {
  xcomplex newz; mpc_sub(newz.z,a.z,b.z,xcomplex::default_rnd); return newz;
}
xcomplex const operator * (const xcomplex &a, const xcomplex &b) {
  xcomplex newz; mpc_mul(newz.z,a.z,b.z,xcomplex::default_rnd); return newz;
}
xcomplex const operator / (const xcomplex &a, const xcomplex &b) {
  xcomplex newz; mpc_div(newz.z,a.z,b.z,xcomplex::default_rnd); return newz;
}

static const xreal xabs(const xcomplex &newz){ 
  xreal tmp;
#ifdef MPFR_REALS
  mpc_abs(tmp.z,newz.z,xcomplex::default_rnd);
#else
  mpfr_t tmpfr;
  mpfr_init2(tmpfr,default_prec);
  mpc_abs(tmpfr,newz.z,xcomplex::default_rnd);
  tmp = mpfr_get_d(tmpfr,xcomplex::default_rnd);
  mpfr_clear(tmpfr);
#endif
  return tmp;
}

static const xreal xabs(const xreal &newz){
#ifdef MPFR_REALS
  xreal tmp;
  mpfr_abs(tmp.z,newz.z,xreal::default_rnd);
  return tmp;
#else
  return fabs(newz);
#endif
}

static const xreal xnorm(const xcomplex &newz){
  xreal tmp;
#ifdef MPFR_REALS
  mpc_norm(tmp.z,newz.z,xcomplex::default_rnd);
#else
  mpfr_t tmpfr;
  mpfr_init2(tmpfr,default_prec);
  tmp = mpfr_get_d(tmpfr,xcomplex::default_rnd);
  mpfr_clear(tmpfr);
#endif
  return tmp;
}

static const long int xlogb(const xcomplex &newz){
  long e = INT_MIN;
  if (mpfr_cmp_si(mpc_realref(newz.z),0) != 0)
    e = mpfr_get_exp(mpc_realref(newz.z));
  if (mpfr_cmp_si(mpc_imagref(newz.z),0) != 0) {
    long ee = mpfr_get_exp(mpc_imagref(newz.z));
    if (ee > e) e = ee;
  }
  return e;
}

static void xscalbln(xcomplex *z, long int a){
  mpfr_mul_2si(mpc_realref(z->z),mpc_realref(z->z),a,xcomplex::default_rnd);
  mpfr_mul_2si(mpc_imagref(z->z),mpc_imagref(z->z),a,xcomplex::default_rnd);
}

static const xreal xroot(const xreal &x, int n)
{
#ifdef MPFR_REALS
  xreal tmp;
#if MPFR_VERSION_MAJOR >= 4
  mpfr_rootn_ui
#else
    mpfr_root
#endif
    (tmp.z,x.z,n,xreal::default_rnd);
  return tmp;
#else
  return pow(x,1.0/n);
#endif
}

static const xreal xsqrt(const xreal &x)
{
#ifdef MPFR_REALS
  xreal tmp;
  mpfr_sqrt(tmp.z,x.z,xreal::default_rnd);
  return tmp;
#else
  return sqrt(x);
#endif
}

static const int xbits(const xcomplex &z)
{
  return default_prec;
  return mpc_get_prec(z.z);
}

static const xreal xeta(const xcomplex &z)
{
#ifdef MPFR_REALS
  return xreal(1,1-default_prec);
#else
  return scalbln(1.0,1-default_prec);
  return scalbln(1.0,1-mpc_get_prec(z.z));
#endif
}

#include "cpoly.C"