1 // arith.h: declarations of arithmetic functions (single precision)
2 //////////////////////////////////////////////////////////////////////////
3 //
4 // Copyright 1990-2012 John Cremona
5 //
6 // This file is part of the eclib package.
7 //
8 // eclib is free software; you can redistribute it and/or modify it
9 // under the terms of the GNU General Public License as published by the
10 // Free Software Foundation; either version 2 of the License, or (at your
11 // option) any later version.
12 //
13 // eclib is distributed in the hope that it will be useful, but WITHOUT
14 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 // for more details.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with eclib; if not, write to the Free Software Foundation,
20 // Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA
21 //
22 //////////////////////////////////////////////////////////////////////////
23
24 #ifndef _ECLIB_ARITH_H
25 #define _ECLIB_ARITH_H 1
26 //flags that this file has been included
27 #include <eclib/interface.h>
28 #include <cstring> // for memset gcc >= 4.3
29 #include <eclib/xmod.h> // supercedes the macros which used to be here
30
31 /* Prime number class; adapted from Pari */
32
33 typedef unsigned char *byteptr;
34
35 class primeclass {
36 friend class primevar;
37 byteptr pdiffptr;
38 long NPRIMES, BIGGESTPRIME;
39
40 void reset(void);
41 int at_end(void);
42 int advance(void);
43 byteptr p_aptr; // points to "current" prime
44 long p_ind; // index of "current" prime
45 long p_val; // value of "current" prime
46
47 public:
48 primeclass(); // will use 10^6 as default or read from file
49 primeclass(long maxnum);
50 ~primeclass();
51 void init(long maxnum); // called in constructor, or to make more primes
52 long number(long n) ; // returns n'th prime (n=1 gives p=2)
53 vector<long> getfirst(long n); // return primes 2..p_n as vector<long>
54 friend long nprimes(void);
55 friend long maxprime(void);
56 };
57
58 extern primeclass the_primes; // The one and only instance
59
prime_number(long n)60 inline long prime_number (long n) /* returns n'th prime from global list */
61 {return the_primes.number(n);}
primes(long n)62 inline vector<long> primes (long n) /* returns list of first n primes */
63 {return the_primes.getfirst(n);}
nprimes(void)64 inline long nprimes(void) {return the_primes.NPRIMES;}
maxprime(void)65 inline long maxprime(void) {return the_primes.BIGGESTPRIME;}
66 long prime_pi(long p); // returns i>=0 such that p is the i'th prime
67
68 class primevar {
69 public:
70 long val; /* current value */
71 long ind; /* current index */
72 private:
73 byteptr ndiff; /* pointer to next diff*/
74 long maxindex; /* max index */
75 public:
76 primevar(long max=the_primes.NPRIMES, long i=1)
77 {maxindex=max; ind=i; val=the_primes.number(i);
78 ndiff=the_primes.pdiffptr+i;}
79 void init(long max=the_primes.NPRIMES, long i=1)
80 {maxindex=max; ind=i; val=the_primes.number(i);
81 ndiff=the_primes.pdiffptr+i;}
82 void operator++() {if ((ind++)<=maxindex) { val+=*ndiff++;}}
83 void operator++(int) {if ((ind++)<=maxindex) { val+=*ndiff++;}}
ok()84 int ok() const {return ind<=maxindex;}
more()85 int more() const {return ind<maxindex;}
value()86 long value() const {return val;}
index()87 long index() const {return ind;}
88 operator long() const {return val;}
89 };
90
91
92 /* Usage of primevar:
93 ---To loop through first n primes:
94
95 long p;
96 for(primevar pr(n); pr.ok(); pr++) {p=pr; ... ;} // or:
97 for(pr.init(n); pr.ok(); pr++) {p=pr; ... ;} // iff pr is existing primevar
98
99 ---To loop through all primes:
100
101 primevar pr; //or for existing primevar, pr.init();
102 long p;
103 while(pr.ok()) {p=pr; pr=++; ...;}
104
105 */
106
107 long primdiv(long); /* returns first prime divisor */
108 vector<long> pdivs(long); /* list of prime divisors */
109
110 vector<long> posdivs(long, const vector<long>& plist); // all positive divisors
posdivs(long n)111 inline vector<long> posdivs(long n)
112 {
113 const vector<long>& plist = pdivs(n);
114 return posdivs(n,plist);
115 }
116
117 vector<long> alldivs(long, const vector<long>& plist); // absolutely all divisors
alldivs(long n)118 inline vector<long> alldivs(long n)
119 {
120 const vector<long>& plist = pdivs(n);
121 return alldivs(n,plist);
122 }
123
124 vector<long> sqdivs(long, const vector<long>& plist); // divisors whose square divides
sqdivs(long n)125 inline vector<long> sqdivs(long n)
126 {
127 const vector<long>& plist = pdivs(n);
128 return sqdivs(n,plist);
129 }
130
131 vector<long> sqfreedivs(long, const vector<long>& plist); // square-free divisors
sqfreedivs(long n)132 inline vector<long> sqfreedivs(long n)
133 {
134 const vector<long>& plist = pdivs(n);
135 return sqfreedivs(n,plist);
136 }
137
138
139 long mod(long a, long modb); /* modulus in range plus or minus half mod */
140 long posmod(long a, long modb); /* ordinary modulus, but definitely positive */
141
142 long gcd(long, long);
143 int gcd(int, int);
144 long lcm(long, long);
145 long bezout(long, long, long&, long&);
146 int intbezout(int aa, int bb, int& xx, int& yy);
147
148 long invmod(long, long);
149 int modrat(long, long, float, long&, long&);
150 int modrat(int, int, float, int&, int&);
151
is_zero(long n)152 inline int is_zero(long n) {return n==0;}
153
154 long val(long factor, long number); // order of factor in number
155
divides(long factor,long number)156 inline int divides(long factor,long number)
157 {
158 return (factor==0) ? (number==0) : (number%factor==0);
159 }
160
ndivides(long factor,long number)161 inline int ndivides(long factor,long number)
162 {
163 return (factor==0) ? (number!=0) : (number%factor!=0);
164 // return !::divides(factor,number);
165 }
166
m1pow(long a)167 inline long m1pow(long a) {return (a%2 ? -1 : +1);}
sign(long a)168 inline int sign(long a) {return (a==0? 0: (a>0? 1: -1));}
sign(double a)169 inline int sign(double a) {return (a==0? 0: (a>0? 1: -1));}
170
171 long chi2(long a);
172 long chi4(long a);
173 long hilbert2(long a, long b);
174 long legendre(long a, long b);
175 long kronecker(long d, long n);
176 int intlog2(long& n, long& e, int roundup);
177
178 int is_squarefree(long n);
179 int is_valid_conductor(long n);
180
181
182 #endif
183