- CLN-list - ginac.de

Floating point overflow discrepancy
by Michael Miller 28 Feb '11

28 Feb '11

The following code (to calculate x^100000): cl_R x="1.0L100000", y=1; for (int i=1; i<=100000; i++) y=x*y; cout << y << "\n"; cout << exp(100000*ln(x)) << "\n";; results in the output 1.00000000000000302104L10000000000 terminate called after throwing an instance of 'cln::floating_point_overflow_exception' what(): floating point overflow. Aborted I know that these are very large numbers, but the two calculations are computing the same value. If the first doesn't overflow, then why does the second? Mike

2 3

Re: [CLN-list] CLN extern "C" patch
by Richard B. Kreckel 15 Feb '11

15 Feb '11

Hi! On 02/15/2011 09:45 AM, Jan wrote: > as far as I know there is no version of MSVC that compiles CLN without > problems to date. I am working on it... Has a bugreport been submitted to the compiler manufacturer regarding this nonconformity with the standard? -richy. -- Richard B. Kreckel <http://www.ginac.de/~kreckel/>

1 0

CLN patch 1
by Jan private 15 Feb '11

15 Feb '11

Problem: MSVC creates different symbol names for 'extern "C" uint32 mulu32_high' when declared inside the cln namespace, and when defined inside it. Linking produces "unresolved external symbol" errors. ----- Which is allowed by the Standard (specifically, [dcl.link]). However, it also says: "At most one function with a particular name can have C language linkage. Two declarations for a function with C language linkage with the same function name (ignoring the namespace names that qualify it) that appear in different namespace scopes refer to the same function. Two declarations for an object with C language linkage with the same name (ignoring the namespace names that qualify it) that appear in different namespace scopes refer to the same object. [Note: because of the one definition rule (3.2), only one definition for a function or object with C linkage may appear in the program; that is, sucha function or object must not be defined in more than one namespace scope." ----- Solution: Moved declarations into global namespace and prefixed with cln_

4 6

mingw-w64
by Jan 08 Feb '11

08 Feb '11

Hi, I tried to compile CLN with MinGW 64-bit (compilation for 32bit works fine). Downloaded mingw-w64-1.0-bin_i686-linux_20110117.tar.bz2 from sourceforge and installed it. git clone git://www.ginac.de/cln.git cd cln autoreconf -i ./configure --host=x86_64-w64-mingw32 --prefix=/usr/x86_64-w64-mingw32 --without-gmp make and this is the result (complete output): Making all in src make[1]: Entering directory `/home/jan/mingw/cln/src' /bin/bash ../libtool --tag=CXX --mode=compile x86_64-w64-mingw32-g++ -DHAVE_CONFIG_H -I. -I../autoconf -I../include -I../src -I../include -I../src -g -O2 -MT cl_alloca.lo -MD -MP -MF .deps/cl_alloca.Tpo -c -o cl_alloca.lo `test -f 'base/cl_alloca.cc' || echo './'`base/cl_alloca.cc libtool: compile: x86_64-w64-mingw32-g++ -DHAVE_CONFIG_H -I. -I../autoconf -I../include -I../src -I../include -I../src -g -O2 -MT cl_alloca.lo -MD -MP -MF .deps/cl_alloca.Tpo -c base/cl_alloca.cc -o cl_alloca.o In file included from ./base/cl_sysdep.h:10:0, from base/cl_alloca.cc:4: ../include/cln/intparam.h:26:2: error: #error "Type char * does not fit into a long!!" In file included from ../include/cln/types.h:10:0, from ./base/cl_macros.h:6, from ./base/cl_sysdep.h:56, from base/cl_alloca.cc:4: ../include/cln/intparam.h:26:2: error: #error "Type char * does not fit into a long!!" base/cl_alloca.cc: In function 'cln::cl_alloca_header* cln::cl_alloc_alloca_header(size_t)': base/cl_alloca.cc:20:40: error: cast from 'long int*' to 'long int' loses precision make[1]: *** [cl_alloca.lo] Error 1 make[1]: Leaving directory `/home/jan/mingw/cln/src' make: *** [all-recursive] Error 1 Any ideas? Thanks, Jan

3 5

[patch] fix compile errors with clang
by Pawel Worach 02 Feb '11

02 Feb '11

Hi, Here are two patches to fix compile errors with clang. libtool: compile: clang++ -DHAVE_CONFIG_H -I. -I../autoconf -I../include -I../src -I../include -I../src -I/usr/local/include -O2 -pipe -fPIC -fno-strict-aliasing -MT cl_rcpointer2_hashweak_rcpointer.lo -MD -MP -MF .deps/cl_rcpointer2_hashweak_rcpointer.Tpo -c base/hash/cl_rcpointer2_hashweak_rcpointer.cc -fPIC -DPIC -o .libs/cl_rcpointer2_hashweak_rcpointer.o In file included from base/hash/cl_rcpointer2_hashweak_rcpointer.cc:7: In file included from ./base/hash/cl_rcpointer2_hashweak_rcpointer.h:8: In file included from ./base/hash/cl_hash2weak.h:6: ./base/hash/cl_hash2.h:121:13: error: use of undeclared identifier '_garcol_fun' if (_garcol_fun(this)) ^ this-> ./base/hash/cl_hash2.h:77:9: note: in instantiation of member function 'cln::cl_heap_hashtable_2<cln::cl_rcpointer, cln::cl_rcpointer, cln::cl_rcpointer>::prepare_store' requested here prepare_store(); ^ base/hash/cl_rcpointer2_hashweak_rcpointer.cc:47:67: note: in instantiation of member function 'cln::cl_heap_hashtable_2<cln::cl_rcpointer, cln::cl_rcpointer, cln::cl_rcpointer>::put' requested here ((cl_heap_weak_hashtable_from_rcpointer2_to_rcpointer*)pointer)->put(x,y,z); ^ In file included from base/hash/cl_rcpointer2_hashweak_rcpointer.cc:7: In file included from ./base/hash/cl_rcpointer2_hashweak_rcpointer.h:8: In file included from ./base/hash/cl_hash2weak.h:6: In file included from ./base/hash/cl_hash2.h:6: ./base/hash/cl_hash.h:36:12: note: must qualify identifier to find this declaration in dependent base class bool (*_garcol_fun) (cl_heap*); // Function to make room in the table. ^ 1 error generated. --- src/base/hash/cl_hash2.h.orig 2009-05-10 22:32:30.000000000 +0200 +++ src/base/hash/cl_hash2.h 2011-01-30 01:59:01.000000000 +0100 @@ -118,7 +118,7 @@ if (this->_freelist < -1) return; // Can we make room? - if (_garcol_fun(this)) + if (this->_garcol_fun(this)) if (this->_freelist < -1) return; // No! Have to grow the hash table. @@ -126,7 +126,7 @@ #else // workaround Sun C++ 4.1 inline function compiler bug if (this->_freelist >= -1) { - if (!_garcol_fun(this) || (this->_freelist >= -1)) + if (!this->_garcol_fun(this) || (this->_freelist >= -1)) grow(); } #endif --- src/base/hash/cl_hashuniq.h.orig 2011-01-30 01:57:23.000000000 +0100 +++ src/base/hash/cl_hashuniq.h 2011-01-30 01:58:31.000000000 +0100 @@ -115,7 +115,7 @@ if (this->_freelist < -1) return; // Can we make room? - if (_garcol_fun(this)) + if (this->_garcol_fun(this)) if (this->_freelist < -1) return; // No! Have to grow the hash table. @@ -123,7 +123,7 @@ #else // workaround Sun C++ 4.1 inline function compiler bug if (this->_freelist >= -1) { - if (!_garcol_fun(this) || (this->_freelist >= -1)) + if (!this->_garcol_fun(this) || (this->_freelist >= -1)) grow(); } #endif -- Pawel

2 1

[cppcheck] Member variable 'gcd_result::gcd_degree' is not initialised in the constructor.
by Martin Ettl 23 Jan '11

23 Jan '11

Hi, during a check with the static code analysis tool cppcheck of the current git head of cln, the tool discoverd a not initialized member variable of class gcd_result::gcd_degree. Please refer the attached patch that fixes the issue. Best regards from the cppcheck team Ettl Martin

2 1

exponants of floating points (or do I need it ?)
by Franck Sabatié 14 Jan '11

14 Jan '11

Hello, I have a hard time seeing how to best implement the following : I have a complicated equation of 4 variables that I need evaluated with a certain number of digits in precision (30 or so) : dummy = (15/4)*(1/(pow(xi,5)*(1+i-k)*(2+i-k)*(3+i-k)*(4+i-k)*(5+i-k)))*( ( (3*pow(3+i-k,2)-(1+i-k)*(5+i-k))*pow(xi*xi-xb,2) -2*(2+i-k)*(4+i-k)*(1-xi*xi)*(xb*xb-xi*xi) ) * ( pow((xb+xi)/(1+xi),3+i-k)- pow((xb-xi)/(1-xi),3+i-k) ) + 6*(3+i-k)*xi*(1-xb)*(xi*xi-xb)*(pow((xb+xi)/(1+xi),3+i-k)+ pow((xb-xi)/(1-xi),3+i-k)) ); where the variables are i, k, xb and xi. Although most of the equation is relatively easy to implement in CLN, I have a hard time figuring out how to deal with the "pow". There appears to be no such exponant function for floats ? Also, it is still unclear to me how to enter this equation painlessly into my code. Do I have to put cl_float(1,30) when I use 1 for instance ? (same for other numbers). Should I evaluate this in floating precision or with "real numbers" ? And finally how do I evaluate expressions such as (something)^(something else) with the (something else) NOT being an integer. Sorry about the very basic questions but most of the demo examples deal with integers and the few that do not, have very little information in them. Thank you for your time, Franck

2 1

Power function
by Joshua Friedman 14 Jan '11

14 Jan '11

I think you need to compute exp(k ln(x) ) for x^k Sent from my iPad On Jan 14, 2011, at 6:00 AM, cln-list-request(a)ginac.de wrote: > Send CLN-list mailing list submissions to > cln-list(a)ginac.de > > To subscribe or unsubscribe via the World Wide Web, visit > https://www.cebix.net/mailman/listinfo/cln-list > or, via email, send a message with subject or body 'help' to > cln-list-request(a)ginac.de > > You can reach the person managing the list at > cln-list-owner(a)ginac.de > > When replying, please edit your Subject line so it is more specific > than "Re: Contents of CLN-list digest..." > > > Today's Topics: > > 1. exponants of floating points (or do I need it ?) (Franck Sabati?) > > > ---------------------------------------------------------------------- > > Message: 1 > Date: Fri, 14 Jan 2011 09:40:23 +0100 > From: Franck Sabati? <fsabatie(a)gmail.com> > To: cln-list(a)ginac.de > Subject: [CLN-list] exponants of floating points (or do I need it ?) > Message-ID: > <AANLkTinG7GBxWw0We%3D8tGv_DNW7CoHK3L-UcFNqaU7iF(a)mail.gmail.com> > Content-Type: text/plain; charset="iso-8859-1" > > Hello, > > I have a hard time seeing how to best implement the following : I have a > complicated equation of 4 variables that I need evaluated with a certain > number of digits in precision (30 or so) : > > dummy = (15/4)*(1/(pow(xi,5)*(1+i-k)*(2+i-k)*(3+i-k)*(4+i-k)*(5+i-k)))*( > ( (3*pow(3+i-k,2)-(1+i-k)*(5+i-k))*pow(xi*xi-xb,2) > -2*(2+i-k)*(4+i-k)*(1-xi*xi)*(xb*xb-xi*xi) ) * ( pow((xb+xi)/(1+xi),3+i-k)- > pow((xb-xi)/(1-xi),3+i-k) ) + > 6*(3+i-k)*xi*(1-xb)*(xi*xi-xb)*(pow((xb+xi)/(1+xi),3+i-k)+ > pow((xb-xi)/(1-xi),3+i-k)) ); > > > where the variables are i, k, xb and xi. Although most of the equation is > relatively easy to implement in CLN, I have a hard time figuring out how to > deal with the "pow". There appears to be no such exponant function for > floats ? Also, it is still unclear to me how to enter this equation > painlessly into my code. Do I have to put cl_float(1,30) when I use 1 for > instance ? (same for other numbers). Should I evaluate this in floating > precision or with "real numbers" ? And finally how do I evaluate expressions > such as (something)^(something else) with the (something else) NOT being an > integer. > > Sorry about the very basic questions but most of the demo examples deal with > integers and the few that do not, have very little information in them. > > Thank you for your time, > Franck >

1 0

Re: [CLN-list] CLN-list Digest, Vol 55, Issue 2
by Joshua Friedman 10 Jan '11

10 Jan '11

Is there an easy to to tell if I succeeded in compiling cln to 64 bit on OS X? On Sat, Dec 11, 2010 at 6:00 AM, <cln-list-request(a)ginac.de> wrote: > Send CLN-list mailing list submissions to > cln-list(a)ginac.de > > To subscribe or unsubscribe via the World Wide Web, visit > https://www.cebix.net/mailman/listinfo/cln-list > or, via email, send a message with subject or body 'help' to > cln-list-request(a)ginac.de > > You can reach the person managing the list at > cln-list-owner(a)ginac.de > > When replying, please edit your Subject line so it is more specific > than "Re: Contents of CLN-list digest..." > > > Today's Topics: > > 1. Re: Does someone see a way to speed this up (Bruno Haible) > > > ---------------------------------------------------------------------- > > Message: 1 > Date: Sat, 11 Dec 2010 00:50:47 +0100 > From: Bruno Haible <bruno(a)clisp.org> > To: cln-list(a)ginac.de > Cc: Joshua Friedman <crowneagle(a)gmail.com> > Subject: Re: [CLN-list] Does someone see a way to speed this up > Message-ID: <201012110050.49085.bruno(a)clisp.org> > Content-Type: text/plain; charset="utf-8" > > Hi Joshua, > > > I am using a lot of function calls, would passing constant pointers speed > up > > my program significantly? ... > > I don't expect people to read the code, just look at the way I am using > cln. > > Is this an efficient way? > > Based on past experience, not on your program in particular, here's a list > of suggestions, starting with the most promising ones, and the least > promising ones at the end. > > 1) Use a profiler that displays the bottlenecks and relative consumption > of CPU time. On Linux/x86, the tool of choice is > valgrind <http://valgrind.org/docs/manual/cl-manual.html>, together > with kcachegrind for the visualization of the results. > > 2) Are you using excess precision in large parts of the computations? > Often, in order to get 20 digits of a result, you need 40 digits > throughout the computation if you choose a constant precision. > But sometimes you can actually do part of the computation with > 22 digits and part with 40 digits, and the result will still have > 20 digits of accuracy. > For analyzing this, you need to have a certain understanding of the > numerical stability of your algorithm. For example, when you are > evaluating integrals along a path around a pole, you use 40 > digits in the computation and get 20 digits accuracy in the result; > but when you are evaluating an integrals along a path around a > region with no singularity, then 40 digits in the computation normally > lead to slightly less than 40 digits accuracy in the result. > > 3) When computing integrals, use higher-order integration methods than > simple averaged summation, if the integrand has continuous derivates. > If your integrand is infinitely differentiable, then Runge-Kutta > (the formula with 1/6, 4/6, 1/6) is a big win, and higher-order > Runge-Kutta methods are even better. > > 4) Pull calls to transcendental functions out of loops if possible. > In CLN, elementary operations and sqrt are fast, whereas sin, exp, > log, etc. are significantly slower. > > 5) Pull common computations out of loops. The compiler won't do it > for it. Things like cl_float(2,prec)*PI are quite fast to compute, > but you may save a couple of percent by moving all these "trivial" > common subexpressions outside the loops. > > 6) Then only, start to worry about whether you pass numbers by > "copy" (e.g. cl_N x) or by reference (e.g. const cl_N& x). > CLN numbers are references anyway, so by doing this you can save > the increment and decrement of a reference count, nothing else. > If you had many function calls and many small-sized bignums > and floats, you could save maybe at most 20% by using this trick. > It's more a micro-optimization than the other five points. > > Bruno > > > ------------------------------ > > _______________________________________________ > CLN-list mailing list > CLN-list(a)ginac.de > https://www.cebix.net/mailman/listinfo/cln-list > > > End of CLN-list Digest, Vol 55, Issue 2 > *************************************** > -- Joshua Friedman PhD CrownEagle(a)gmail.com http://www.math.sunysb.edu/~joshua

2 1

Does someone see a way to speed this up
by Joshua Friedman 11 Dec '10

11 Dec '10

I am using a lot of function calls, would passing constant pointers speed up my program significantly? I am doing computations with modular forms (2 variable functions, written as a Fourier series in one variable). I am basically computing integrals of these functions and searching for maximums. The program takes a long time to run, I think much longer than it should. I don't expect people to read the code, just look at the way I am using cln. Is this an efficient way? #include <iostream.h> #include <fstream.h> #include <stdio.h> #include <stdlib.h> #include <string.h> #include <cln/real.h> #include <cln/integer.h> #include <iostream.h> #include <fstream.h> #include <cln/io.h> #include <cln/integer_io.h> #include <cln/real_io.h> #include <cln/float_io.h> #include <cln/float_io.h> #include <cln/integer.h> #include <cln/real.h> #include <cln/complex.h> // We do I/O. #include <cln/io.h> #include <cln/integer_io.h> #include <cln/real_io.h> #include <cln/complex_io.h> // We use the timing functions. #include <cln/timing.h> //using namespace std; using namespace cln; typedef cln::cl_I bigI; typedef cln::cl_F bigF; typedef cln::cl_N bigC; typedef cln::cl_R bigR; bigF h(bigF y, bigF a[], int k, int M, int pr); bigF ff(bigF x, bigF y, bigF a[], int k, int M, int pr); bigF v(bigF y, bigF a[], int k, int M, int pr); bigF integralH(bigF a[], int k, int M, int n, int pr); bigF integralL(bigF a[], int k, int M, int n, int pr); bigF max(bigF a[][100],bigF b[], int k, int D, int M, bigF e, int pr); //for each weight k (2k gives actual weight), make 2d arrary of coefficients //write routine to search FD domain (less than 2k/pi, check the 2), first do //sup for indiv forms, than the average one. int main(int argc, char *argv[]) { //the number of coefficients to use int pbits= 128; int nPartitions = 100; int divisor = 4; float fbinsize = 0.1; int startK = 22; cout << "\n Enter prec, numPartitions, intervalSize, divisor to reduce coef, startK "; cin >> pbits >> nPartitions >> fbinsize >> divisor >> startK; std::string snum; char w; int K,P,D; bigI S, H; std::ifstream myfile ("full_output.txt"); if (myfile.is_open()) { while (! myfile.eof() ) { D = 0; myfile >> w >> K >> P >> D; bigI bigInt; int pr = pbits; cln::float_format_t prec = cln::float_format_t(pr); bigF binSize = cl_float(fbinsize,prec); bigF a[D][100]; bigF b[D]; { for (int n = 0; n < D; n++) { for (int u = 0; u < P; u++) { myfile >> S >> bigInt; a[n][u] = cl_float(bigInt,prec); } if (K >= startK) { b[n] = integralH(a[n],K, P/divisor, nPartitions, pr) + integralL(a[n],K, P/divisor, nPartitions, pr); } //b[n] = cl_float(1,prec); //cout << "\nLoaded Coeff"; } if (K >= startK) { if (D>0) { bigF sNorm = max(a,b,K,D,P, binSize, pr); cout << "\n" << K << "," << sNorm; } } } } myfile.close(); } else cout << "Unable to open file"; return 0; } //Used to do the integral over cylindar y > 1 bigF h(bigF y, bigF a[], int k, int M, int pr) { cln::float_format_t prec = cln::float_format_t(pr); const bigF PI = pi(prec); bigF s = cln::cl_float(0,prec); for (int n = 1; n < M; n++) if (cl_float(4,prec)*PI*cl_float(n,prec)*y < cl_float(300,prec)) s = s + a[n]*a[n]*exp(cl_float(-4,prec)*PI*cl_float(n,prec)*y); return exp(ln(y)*cl_float(2*k-2, prec))*s; } //ff is equiv to g, use it as a test and see which is more effiecient bigF ff(bigF x, bigF y, bigF a[], int k, int M, int pr) { cln::float_format_t prec = cln::float_format_t(pr); const bigF PI = pi(prec); bigF q = cl_float(0.0,prec); int j = 1; bigF w = cl_float(0.0,prec); int r = 0; bigF hh = cl_float(0,prec); for (j = 1; j <= M; j++) { if (cl_float(4,prec)*PI*cl_float(j,prec)*y < cl_float(50,prec)) w = w + a[j]*a[j]*exp(-cl_float(4,prec)*PI*cl_float(j,prec)*y); } for (r = 3; r <= M; r++) { q = cl_float(0.0,prec); for (j = 1; 2*j < r; j++) q = q + a[j]*a[r-j]*cos(cl_float(2,prec)*PI*(cl_float(2,prec)*cl_float(j,prec) - cl_float(r,prec))*x); if (cl_float(2,prec)*PI*cl_float(r,prec)*y < cl_float(50,prec)) q = q*exp(cl_float(-2,prec)*PI*cl_float(r,prec)*y); else q = cl_float(0,prec); w = w + 2*q; } bigF tt = cln::cl_float(2*k,prec); return exp(ln(y)*tt)*w; } bigF max(bigF a[][100], bigF b[], int k, int D, int M, bigF e, int pr) { bigF x,y,half; cln::float_format_t prec = cln::float_format_t(pr); bigF one = cl_float(1,prec); x = cl_float(0,prec); y = cl_float(1,prec); half = cl_float(0.5,prec); bigF K = cl_float(k,prec); bigF mxT = cl_float(0,prec); bigF sum = cl_float(0,prec); bigF w = sqrt(cl_float(3,prec))/cl_float(2,prec); while ( x < half) { y = cl_float(1,prec); while (y < K/cl_float(6,prec)) { sum = cl_float(0,prec); for (int i = 0; i < D; i++) sum = sum + ff(x,y,a[i],k,M,pr)/b[i]; if (sum > mxT) mxT = sum; y = y + e; } x = x + e; } bigF mxB = cl_float(0,prec); x = cl_float(0,prec); while ( x < half) { y = sqrt(cl_float(1) - x*x); while (y < one) { sum = cl_float(0,prec); for (int i = 0; i < D; i++) sum = sum + ff(x,y,a[i],k,M,pr)/b[i]; if (sum > mxB) mxB = sum; y = y + e; } x = x + e; } if (mxB > mxT) return mxB; return mxT; } //twice the Integral f(z) square in x direction from 0 to 1/2; //y is variable from sqrt(3)/2 to 1 bigF v(bigF y, bigF a[], int k, int M, int pr) { cln::float_format_t prec = cln::float_format(pr); const bigF PI = pi(prec); bigF q = cl_float(0.0, prec); bigF w = cl_float(0.0, prec); bigF one = cl_float(1, prec); bigF two = cl_float(2, prec); bigF three = cl_float(3, prec); bigF four = cl_float(4, prec); for (int j = 1; j <= M; j++) { bigF J = cl_float(j,prec); if (four*PI*J*y < cl_float(50,prec)) w = w + a[j]*a[j]*exp(-four*PI*J*y); } w = two*(cl_float(0.5,prec) - sqrt(one-y*y))*w; for (int r = 3; r <= M; r++) { bigF R = cl_float(r,prec); q = cl_float(0,prec); for (int j = 1; 2*j < r; j++) { bigF J = cl_float(j,prec); q = q + a[j]*a[r-j]/(two*PI*(two*J - R))* (sin(PI*(two*J - R)) - sin(two*PI*(two*J - R)*sqrt(one-y*y))); } //q = q + a[j]*a[r-j]*cos(2*PI*(2*j - r)*x); if (two*PI*R*y < cl_float(50,prec)) q = q*exp(-two*PI*R*y); else q = cl_float(0,prec); w = w + two*q; } //for (int n = 1; n < h; n++) s = s + a[n]*exp(2*PI*I*z*n); bigF tt = cln::cl_float((float)(2*k-2),prec); return exp(ln(y)*tt)*w; } bigF integralH(bigF a[], int k, int M, int n, int pr) { const bigF PI = "3.14159265358979323846264338327950288"; cln::float_format_t prec = cln::float_format_t(pr); bigF q = cl_float(0.0, prec); bigF w = cl_float(0.0, prec); bigF aa = cl_float(0.0, prec); bigF b = cl_float(0.0, prec); bigF d = cl_float(0.0, prec); bigF one = cl_float(1, prec); bigF two = cl_float(2, prec); bigF three = cl_float(3, prec); bigF four = cl_float(4, prec); bigF five = cl_float(5, prec); aa = one; //modify this later b = cl_float(k,prec)/two; d = (b-aa)/cl_float(n,prec); //f(x) is the dummy function q = h(aa,a,k,M, pr) + h(b,a,k,M, pr); for (int j = 1; 2*j <= n-2; j++) w = w + h(aa + (two*cl_float(j,prec))*d,a,k,M, pr); q = q + two*w; w = cl_float(0,prec); for (int j = 1; 2*j <= n; j++) w = w + h(aa + (two*cl_float(j,prec)-one)*d,a,k,M, pr); q = q + four*w; q = q*d/three; return q; } bigF integralL(bigF a[],int k, int M, int n, int pr) { cln::float_format_t prec = cln::float_format(pr); bigF q = cl_float(0.0, prec); bigF w = cl_float(0.0, prec); bigF aa = cl_float(0.0, prec); bigF b = cl_float(1.0, prec); bigF d = cl_float(0.0, prec); bigF one = cl_float(1, prec); bigF two = cl_float(2, prec); bigF three = cl_float(3, prec); bigF four = cl_float(4, prec); bigF five = cl_float(5, prec); aa = sqrt(three)/two; b = one; d = (b-aa)/cl_float(n,prec); //f(x) is the dummy function q = v(aa,a,k,M, pr) + v(b,a,k,M, pr); for (int j = 1; 2*j <= n-2; j++) w = w + v(aa + two*cl_float(j,prec)*d,a,k,M, pr); q = q + two*w; w = cl_float(0,prec); for (int j = 1; 2*j <= n; j++) w = w + v(aa + (two*cl_float(j,prec)-one)*d,a,k,M, pr); q = q + four*w; q = q*d/three; return q; } -- Joshua Friedman PhD CrownEagle(a)gmail.com http://www.math.sunysb.edu/~joshua

2 1