← Back to branch summary

~ubuntu-branches/ubuntu/quantal/gnumeric/quantal

« back to all changes in this revision

Viewing changes to plugins/numtheory/numtheory.c

  • Committer: Bazaar Package Importer
  • Author(s): Gauvain Pocentek
  • Date: 2009-06-07 11:10:47 UTC
  • mfrom: (1.1.19 upstream) (2.1.2 squeeze)
  • Revision ID: james.westby@ubuntu.com-20090607111047-l3rtbzfjxvmi1kx0
Tags: 1.9.8-1ubuntu1
* Merge from debian unstable, remaining changes:
  - Promoted gnumeric-doc to Recommends in gnumeric package for help to be
    installed automatically
  - gnumeric-gtk is a transitional package
  - gnumeric conflicts with gnumeric-gtk << 1.8.3-3ubuntu1
  - call initltool-update in po*
  - remove psiconv support (psiconv is in universe):
    o debian/control: remove B-D on libpsiconv-dev
    o debian/rules: don't pass --with-psiconv to ./configure

Show diffs side-by-side

added added

removed removed

Lines of Context:
plugins/numtheory/numtheory.c
33
33
 
34
34
#define OUT_OF_BOUNDS "#LIMIT!"
35
35
 
 
36
/*
 
37
 * The largest integer i, such at all integers {0,...,i} can be accurately
 
38
 * represented in a gnm_float _and_ in a guint64.  (For regular "double",
 
39
 * the latter part is irrelevant.)
 
40
 */
36
41
static const double bit_max = MIN (1 / GNM_EPSILON, (gnm_float)G_MAXUINT64);
37
42
 
38
43
/* ------------------------------------------------------------------------- */
39
44
 
40
 
static int
 
45
static GnmValue *
 
46
value_new_guint64 (guint64 n)
 
47
{
 
48
        return value_new_float (n);
 
49
}
 
50
 
 
51
 
 
52
static guint64
41
53
intpow (int p, int v)
42
54
{
43
 
        int temp;
 
55
        guint64 temp;
44
56
 
45
57
        if (v == 0) return 1;
46
58
        if (v == 1) return p;
50
62
        return (v % 2) ? temp * p : temp;
51
63
}
52
64
 
53
 
#define PTABLE_CHUNK 64
54
 
#define ITHPRIME_LIMIT (1 << 20)
 
65
#define ITHPRIME_LIMIT (1 << 22)
55
66
static gint *prime_table = NULL;
56
67
 
57
68
/* Calculate the i-th prime.  Returns TRUE on error.  */
58
69
static gboolean
59
 
ithprime (int i, int *res)
 
70
ithprime (int i, guint64 *res)
60
71
{
61
72
        static int computed = 0;
62
73
        static int allocated = 0;
68
79
                int candidate;
69
80
 
70
81
                if (i > allocated) {
71
 
                        g_assert (PTABLE_CHUNK >= 2);
72
 
                        allocated = MAX (i, allocated + PTABLE_CHUNK);
 
82
                        allocated = MAX (i, 2 * allocated + 100);
 
83
                        allocated = MIN (allocated, ITHPRIME_LIMIT);
73
84
                        prime_table = g_renew (int, prime_table, allocated);
74
85
                        if (computed == 0) {
75
86
                                prime_table[computed++] = 2;
107
118
 * Returns TRUE on error.
108
119
 */
109
120
static gboolean
110
 
walk_factorization (int n, void *data,
111
 
                    void (*walk_term) (int p, int v, void *data))
 
121
walk_factorization (guint64 n, void *data,
 
122
                    void (*walk_term) (guint64 p, int v, void *data))
112
123
{
113
 
        int index = 1, p = 2, v;
 
124
        int index = 1, v;
 
125
        guint64 p = 2;
114
126
 
115
127
        while (n > 1 && p * p <= n) {
116
128
                if (ithprime (index, &p))
144
156
/*
145
157
 * Returns -1 (out of bounds), or #primes <= n
146
158
 */
147
 
static int
148
 
compute_nt_pi (int n)
 
159
static gint64
 
160
compute_nt_pi (guint64 n)
149
161
{
150
 
        int lower = 2, upper = 4, mid, p = 7;
 
162
        guint64 lower = 2, upper = 4, mid, p = 7;
151
163
 
152
164
        if (n <= 1)
153
165
                return 0;
184
196
static int
185
197
isprime (guint64 n)
186
198
{
187
 
        int i = 1, p = 2;
 
199
        int i = 1;
 
200
        guint64 p = 2;
188
201
 
189
202
        if (n <= 1)
190
203
                return 0;
191
204
 
192
 
        for (i = 1; (guint64)p * (guint64)p <= n; i++) {
 
205
        for (i = 1; p * p <= n; i++) {
193
206
                if (ithprime (i, &p))
194
207
                        return -1;
195
208
                if (n % p == 0)
216
229
};
217
230
 
218
231
static void
219
 
walk_for_phi (int p, int v, void *data)
 
232
walk_for_phi (guint64 p, int v, void *data_)
220
233
{
221
 
        *((int *)data) *= intpow (p, v - 1) * (p - 1);
 
234
        guint64 *data = data_;
 
235
        *data *= intpow (p, v - 1) * (p - 1);
222
236
}
223
237
 
224
238
static GnmValue *
225
239
gnumeric_phi (GnmFuncEvalInfo *ei, GnmValue const * const *args)
226
240
{
227
 
        int phi = 1;
 
241
        guint64 phi = 1;
228
242
        gnm_float n = gnm_floor (value_get_as_float (args[0]));
229
243
 
230
 
        if (n < 1 || n > INT_MAX)
 
244
        if (n < 1 || n > bit_max)
231
245
                return value_new_error_NUM (ei->pos);
232
246
 
233
 
        if (walk_factorization ((int)n, &phi, walk_for_phi))
 
247
        if (walk_factorization ((guint64)n, &phi, walk_for_phi))
234
248
                return value_new_error (ei->pos, OUT_OF_BOUNDS);
235
249
 
236
 
        return value_new_int (phi);
 
250
        return value_new_guint64 (phi);
237
251
}
238
252
 
239
253
/* ------------------------------------------------------------------------- */
258
272
};
259
273
 
260
274
static void
261
 
walk_for_mu (int p, int v, void *data)
 
275
walk_for_mu (guint64 p, int v, void *data_)
262
276
{
263
 
        *((int *)data) = (v >= 2) ?  0 : - *((int *)data) ;
 
277
        int *data = data_;
 
278
        *data = (v >= 2) ?  0 : -*data;
264
279
}
265
280
 
266
281
static GnmValue *
269
284
        int mu = 1;
270
285
        gnm_float n = gnm_floor (value_get_as_float (args[0]));
271
286
 
272
 
        if (n < 1 || n > INT_MAX)
 
287
        if (n < 1 || n > bit_max)
273
288
                return value_new_error_NUM (ei->pos);
274
289
 
275
 
        if (walk_factorization ((int)n, &mu, walk_for_mu))
 
290
        if (walk_factorization ((guint64)n, &mu, walk_for_mu))
276
291
                return value_new_error (ei->pos, OUT_OF_BOUNDS);
277
292
 
278
293
        return value_new_int (mu);
295
310
};
296
311
 
297
312
static void
298
 
walk_for_d (int p, int v, void *data)
 
313
walk_for_d (guint64 p, int v, void *data_)
299
314
{
300
 
        * (int *) data *= (v + 1);
 
315
        int *data = data_;
 
316
        *data *= (v + 1);
301
317
}
302
318
 
303
319
static GnmValue *
306
322
        int d = 1;
307
323
        gnm_float n = gnm_floor (value_get_as_float (args[0]));
308
324
 
309
 
        if (n < 1 || n > INT_MAX)
 
325
        if (n < 1 || n > bit_max)
310
326
                return value_new_error_NUM (ei->pos);
311
327
 
312
 
        if (walk_factorization ((int)n, &d, walk_for_d))
 
328
        if (walk_factorization ((guint64)n, &d, walk_for_d))
313
329
                return value_new_error (ei->pos, OUT_OF_BOUNDS);
314
330
 
315
331
        return value_new_int (d);
331
347
};
332
348
 
333
349
static void
334
 
walk_for_sigma (int p, int v, void *data)
 
350
walk_for_sigma (guint64 p, int v, void *data_)
335
351
{
336
 
        * (int *) data *=
337
 
                  ( v == 1 ? p + 1 : (intpow (p, v + 1) - 1) / (p - 1) );
 
352
        guint64 *data = data_;
 
353
        *data *= ( v == 1 ? p + 1 : (intpow (p, v + 1) - 1) / (p - 1) );
338
354
}
339
355
 
340
356
static GnmValue *
341
357
gnumeric_sigma (GnmFuncEvalInfo *ei, GnmValue const * const *args)
342
358
{
343
 
        int sigma = 1;
 
359
        guint64 sigma = 1;
344
360
        gnm_float n = gnm_floor (value_get_as_float (args[0]));
345
361
 
346
 
        if (n < 1 || n > INT_MAX)
 
362
        if (n < 1 || n > bit_max)
347
363
                return value_new_error_NUM (ei->pos);
348
364
 
349
 
        if (walk_factorization ((int)n, &sigma, walk_for_sigma))
 
365
        if (walk_factorization ((guint64)n, &sigma, walk_for_sigma))
350
366
                return value_new_error (ei->pos, OUT_OF_BOUNDS);
351
367
 
352
 
        return value_new_int (sigma);
 
368
        return value_new_guint64 (sigma);
353
369
}
354
370
 
355
371
/* ------------------------------------------------------------------------- */
370
386
static GnmValue *
371
387
gnumeric_ithprime (GnmFuncEvalInfo *ei, GnmValue const * const *args)
372
388
{
373
 
        int p;
 
389
        guint64 p;
374
390
        gnm_float i = gnm_floor (value_get_as_float (args[0]));
375
391
 
376
392
        if (i < 1 || i > INT_MAX)
379
395
        if (ithprime ((int)i, &p))
380
396
                return value_new_error (ei->pos, OUT_OF_BOUNDS);
381
397
 
382
 
        return value_new_int (p);
 
398
        return value_new_guint64 (p);
383
399
}
384
400
 
385
401
/* ------------------------------------------------------------------------- */
420
436
 
421
437
/*
422
438
 * Returns
423
 
 *   -1 (out of bounds)
424
 
 *    0 (n <= 1)
 
439
 *    0 (n <= 1) or (out of bounds)
425
440
 *    smallest prime facter
426
441
 */
427
 
static int
 
442
static guint64
428
443
prime_factor (guint64 n)
429
444
{
430
 
        int i = 1, p = 2;
 
445
        int i = 1;
 
446
        guint64 p = 2;
431
447
 
432
448
        if (n <= 1)
433
449
                return 0;
434
450
 
435
 
        for (i = 1; (guint64)p * (guint64)p <= n; i++) {
 
451
        for (i = 1; p * p <= n; i++) {
436
452
                if (ithprime (i, &p))
437
 
                        return -1;
 
453
                        return 0;
438
454
                if (n % p == 0)
439
455
                        return p;
 
456
 
440
457
        }
441
458
 
442
459
        return n;
461
478
gnumeric_pfactor (GnmFuncEvalInfo *ei, GnmValue const * const *args)
462
479
{
463
480
        gnm_float n = gnm_floor (value_get_as_float (args[0]));
464
 
        int p;
 
481
        gint64 p;
465
482
 
466
483
        if (n < 2)
467
484
                return value_new_error_VALUE (ei->pos);
468
485
        if (n > bit_max)
469
 
                p = -1;
 
486
                p = 0;
470
487
        else
471
488
                p = prime_factor ((guint64)n);
472
489
 
473
 
        if (p < 0)
 
490
        if (p == 0)
474
491
                return value_new_error (ei->pos, OUT_OF_BOUNDS);
475
492
 
476
 
        return value_new_int (p);
 
493
        return value_new_float (p);
477
494
}
478
495
 
479
496
/* ------------------------------------------------------------------------- */
496
513
gnumeric_nt_pi (GnmFuncEvalInfo *ei, GnmValue const * const *args)
497
514
{
498
515
        gnm_float n = gnm_floor (value_get_as_float (args[0]));
499
 
        int pi;
 
516
        gint64 pi;
500
517
 
501
518
        if (n < 0)
502
519
                pi = 0;
503
 
        else if (n > INT_MAX)
 
520
        else if (n > bit_max)
504
521
                pi = -1;
505
522
        else
506
 
                pi = compute_nt_pi (n);
 
523
                pi = compute_nt_pi ((guint64)n);
507
524
 
508
525
        if (pi == -1)
509
526
                return value_new_error (ei->pos, OUT_OF_BOUNDS);