Division gets funny?

Zhylian (28)
Hi guys, here I wrote a program to evaluate the integral of a given function. Something goes wrong when evaluating the parameter m_h of the class Integrator.hxx when the method 'Integral_Eval ( double a, double b, unsigned int N )' is called.
I checked the value of all the variables involved whit GDB, everything is fine until the command at line 43 of the file MiddlePoint.cxx get hit.
It should be Pi/10 ( 0.314159...) but it is 0,2999...
The funny things is that this problem occours only when a = 0 and b = Pi ( or viceversa), but as a = 0, b = 1 and a = 1, b=2 everything works fine ( see the assert at the beginning of the main.cxx file )

I will leave all the code here just in case someone wants to run the whole thing. Thanks in advance


***main.cxx***
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
 
  #include "Function.hxx"    		//it already include BaseFunction.hxx
#include "MiddlePoint.hxx" 		//it already include Integrator.hxx

#include<iostream>
#include<iomanip>
#include<assert.h>
#include<cmath>

int main () {

	//assert
	Sine sine_test ( 1., 0, 0, 1. ) ;
	MiddlePoint Test ( &sine_test ) ;

	//m_h class member is fine...
	assert ( std::fabs( Test.Integral_Eval ( 0, 1 , 10 ) - 0.45988929071851814 ) < 1.e-7 ) ;
	assert ( std::fabs( Test.Integral_Eval ( 1, 2 , 30 ) - 0.9564934239032155 ) < 1.e-7 ) ;

	//m_h class member is wrong!! WHY?!
	std::cout << Test.Integral_Eval ( M_PI, 0 , 10 ) << std::endl << Test.getH() << std::endl;	
	//assert ( std::fabs( Test.Integral_Eval ( M_PI, 0 , 10 ) + 2.0082484079079745 ) < 1.e-7 ) ;
	std::cout << Test.Integral_Eval ( 0, M_PI , 100 ) << std::endl << Test.getH() << std::endl;
	//assert ( std::fabs( Test.Integral_Eval ( 0, M_PI , 100 ) + 2.000082249070986 ) < 1.e-7 ) ;
																					//1
	std::cout << "evaluation of the integral of sin(x) between 0, pi\n" ;

	double a {} ;
	double b { M_PI } ;

	Sine sine ( 1., 0, 0, 1. ) ;  //sin(x)
	BaseFunction* f = &sine ;
	
	MiddlePoint M ( a, b, 10 ) ;

	double integral = M.Integral_Eval ( f ) ;
	
	std::cout << std::fixed ;
	std::cout << std::setprecision( 10 ) << integral << std::endl ;
																					//2
	a = 0 ;
	b = 1 ;

	MiddlePoint A ( a, b, 10 ) ;
	
	integral = A.Integral_Eval ( f ) ;
	assert ( std::fabs( integral - 0.45988929071851814 ) < 1.e-7 ) ;

	std::cout << "evaluation of the integral of sin(x) between 0, 1\n" ;
	std::cout << integral << std::endl ;

																					//3
	a = 1 ;
	b = 2 ;

	MiddlePoint B ( a, b, 30 ) ;
	
	integral = B.Integral_Eval ( &sine ) ;
	assert ( std::fabs( integral - 0.9564934239032155 ) < 1.e-7 ) ;

	std::cout << "evaluation of the integral of sin(x) between 0, 1\n" ;
	std::cout << integral << std::endl ;
																					//4

	a = M_PI ;
	b = 0 ;

	MiddlePoint C ( a, b, 10 ) ;
	
	integral = C.Integral_Eval ( f ) ;

	std::cout << "evaluation of the integral of sin(x) between pi, 0\n" ;
	std::cout << integral << std::endl ;

}


***Integrator.hxx***
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
 
#pragma once

#include "BaseFunction.hxx"

#include<cmath>



class Integrator {

 public:

 virtual double Integral_Eval ( const BaseFunction* ) = 0 ;  
 virtual double Integral_Eval ( double ,double , unsigned int ) = 0 ;

 											//constructors
 Integrator () = default ;
 
 Integrator ( double a, double b, unsigned int N )  { 
 //no need to do anything if a < b
 if ( a > b ) {
 	m_a = b ;
 	m_b = a ;
 	m_sign = -1 ;
 	} else {
 		m_a = a ;
 		m_b = b ;
 	}
 m_steps = N ;
 m_h = (double) ( m_b - m_a ) / m_steps ;
 } ;
 
 Integrator ( double a, double b ) {
  //no need to do anything if a < b
 	 if ( a > b ) {
 	 	m_a = b ;
 	 	m_b = a ;
 	 	m_sign = -1 ;
 	 	} else {
 	 		m_a = a ;
 	 		m_b = b ;
 	 	}
 	 m_h = ( m_b - m_a ) / m_steps ;	
 };

 Integrator ( BaseFunction* F ) : m_a {} , m_b {} , m_h {} , m_f { F } { } ;


											//virtual destructor
 virtual ~Integrator() = default ;

											//getter
 [[nodiscard]]
 int getA () const { return m_a ; } ;

 [[nodiscard]]
 int getB () const { return m_b ; } ;

 [[nodiscard]]
 int getSign () const { return m_sign ; } ;
 
 [[nodiscard]]
 double getH () const { return m_h ; } ;
 
 [[nodiscard]]
 unsigned int getN_Steps () const { return m_steps ; } ;

 [[nodiscard]]
 BaseFunction* getFunction () const { return m_f ; } ;
 
											//setter
 void setA ( double a ) { m_a = a ; } ;
 void setB ( double b ) { m_b = b ; } ;
 void setSign ( int s ) { m_sign = s ; } ;
 void setN_Steps ( unsigned int N ) { m_steps = N ; } ;
 void setH ( double h ) { m_h = h ; } ;
 void setFunction ( BaseFunction* F ) { m_f = F ; } ;
 
 private:

 double m_a , m_b ;    				
 int m_sign = 1 ;							//default value
 unsigned int m_steps = 100 ;				//default value
 double m_h ;
 BaseFunction* m_f ;
 
} ;


***MiddlePoint.hxx***
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
 
  agma once

#include "Integrator.hxx"
#include "BaseFunction.hxx"

class MiddlePoint : public Integrator {  //private members of Integrator will be get access by getter/setter methods ( of Integrator )

 public: 

 MiddlePoint () = default ;
 //It's a shadow-usage of the Integrator construcor
 MiddlePoint ( double, double, unsigned int ) ; 
 MiddlePoint ( double , double ) ;
 MiddlePoint ( BaseFunction* ) ;

 ~MiddlePoint () = default ;

 [[nodiscard]]
 virtual double Integral_Eval ( const BaseFunction* ) override ;

 [[nodiscard]]
 virtual double Integral_Eval ( double ,double , unsigned int ) override ;

 private:

 double m_sum = 0 ;
 double m_integral ;
	
} ;


***MiddlePoint.cxx***
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
 
 #include "MiddlePoint.hxx" 


MiddlePoint::MiddlePoint ( double a, double b, unsigned int N ) : Integrator ( a, b, N ) { } 

MiddlePoint::MiddlePoint ( double a, double b ) : Integrator ( a, b ) { } 

MiddlePoint::MiddlePoint ( BaseFunction* F ) : Integrator ( F ) { }

double MiddlePoint::Integral_Eval ( const BaseFunction* F ) {
	//initial check
	if ( getA() == getB() ) return 0 ;

	//algoritm
	for ( unsigned int i {} ; i < getN_Steps () ; ++i ) {
		m_sum += F->Evaluate ( getA() + ( i + 0.5 ) * getH() ) ;
		
	}
	m_integral = getSign() * m_sum * getH() ;
	m_sum = 0; //reset to default value

	return m_integral ;
	
}

double MiddlePoint::Integral_Eval ( double a, double b, unsigned int N ) {
	//initial check
	if ( a == b ) return 0 ;

	//build the interval
	if ( a < b ) {
		setA ( a ) ;
		setB ( b ) ;
		setSign ( 1 ) ;
	} else {
		setA ( b ) ;
		setB ( a ) ;
		setSign ( -1 ) ;
	} 

	//set parameters
	setN_Steps ( N ) ;
	setH ( (double) ( getB() - getA() ) / getN_Steps() ) ;
	
	//algoritm
	m_sum = 0 ;
	for ( unsigned int i {} ; i < getN_Steps() ; ++i ) {
		m_sum += getFunction()->Evaluate ( getA() + ( i + 0.5 ) * getH() ) ;
		
	}
	m_integral = getSign() * m_sum * getH() ;
	m_sum = 0; //reset to default value

	return m_integral ;
	
}


***BaseFunction.hxx***
1
2
3
4
5
6
7
8
9
10
11
 
 #pragma once


class BaseFunction {

 public:
 
 virtual double Evaluate ( double x ) const = 0 ;
 virtual ~BaseFunction () = default ;  
 
} ;


***Function.hxx***
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
 
  #pragma once

#include "BaseFunction.hxx" 
#include<cmath>


/****************************************************** SINE */
class Sine : public BaseFunction {

 public:

 Sine () = default ;
 Sine ( double a, double b, double c, double d ) : m_a { a } , m_b { b } , m_c { c } , m_d { d } { } ; //Sine (1,0,0,1) for standard f(x) = sin(x)
 ~Sine () = default ;

 double Evaluate ( double x ) const override { return m_d * sin ( m_a * x + m_b ) + m_c ; } ;

 [[nodiscard]] 
 double getA () const { return m_a ; } ;
 [[nodiscard]]
 double getB () const { return m_b ; } ;
 [[nodiscard]]
 double getC () const { return m_c ; } ;
 [[nodiscard]]
 double getD () const { return m_d ; } ;

 void setA ( double a ) { m_a = a ; } ;
 void setB ( double b ) { m_b = b ; } ;
 void setC ( double c ) { m_c = c ; } ;
 void setD ( double d ) { m_d = d; } ;

 private:

 double m_a , m_b , m_c, m_d ;
 
} ;


***Makefile***
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
 
  CXXFLAGS = -Wall --pedantic -g -std=c++11 
INCS = `root-config --libs`
LIBS = `root-config --cflags`

esercizio8.0: main.o MiddlePoint.o

	g++ -o $@ main.o MiddlePoint.o ${INCS}

main.o: main.cxx Function.hxx
	g++ -c main.cxx ${CXXFLAGS} ${LIBS}

MiddlePoint.o: MiddlePoint.cxx MiddlePoint.hxx Integrator.hxx BaseFunction.hxx
	g++ -c MiddlePoint.cxx ${CXXFLAGS}


clean:
	rm *.o

run:
	./esercizio8.0

open:
	micro main.cxx Integrator.hxx MiddlePoint.hxx MiddlePoint.cxx 

debug:
	gdb ./esercizio8.0
Last edited on
salem c (3700)
Few people will bother to copy/paste/repair 500 lines of code spread across half a dozen files.
http://sscce.org/

It's up to you to
a) copy all that code into a single source file.
b) start chopping out all the code which isn't necessary to show the problem.

If in the process of removing fluff, you discover the problem - job done.

If you don't find the problem, you'll have a nice <50 line single file demonstration of the issue that people will be far more inclined to study.
TheIdeasMan (6817)
setH ( (double) ( getB() - getA() ) / getN_Steps() ) ;

Hi

The above code at first glance looks like it contains integer division, but that is not the case: A, B, N_Steps are ints. How about coding those get functions to return doubles instead of ints? So what is happening ? When A is M_PI, getA() returns 3 hence the problem.

Another thing that is pedantic, but I find helpful: Always write literal double values with digits either side of the decimal point, as in 10.0 not 10 . In my mind it reinforces the idea that one is dealing with a double; takes away confusion that it looks like an int.

I would recommend the C++ static_cast over the C style cast, it's easier to see and less error prone because the C style cast does no checks - it blindly casts the value, come what may.

With the compilation, have the flag -Wextra as well.

Is there a reason that you are compiling to c++11 standard? c++20 is the current standard.

I am not sure on the details of makefile these days, (my IDE writes it for me, I do edit a couple of things in the cmake file) , but does one have to compile the header files as well? I could compile the whole thing one 1 line from the shell, by compiling only the *.cxx files. Ah I see why now. you have implementation code in your header files - really try to avoid doing that. You could mark them inline but I would personally have implementation code in the *.cxx files, only declarations in the header files.
Last edited on
lastchance (6980)
@Zhylian,
How is it possible that you can make numerical integration so staggeringly complex?

Try this:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 
#include <iostream>
#include <cmath>
using namespace std;

const double PI = 4.0 * atan( 1.0 );

double Midpoint( double f( double ), double a, double b, int N )
{
   double sum = 0;
   double dx = ( b - a ) / N;
   for ( int i = 0; i < N; i++ ) sum += f( a + ( i + 0.5 ) * dx );
   return dx * sum;
}

int main()
{
   int N = 1000;
   double a = 0.0, b = PI;
   auto f = []( double x ){ return sin( x ); };
   cout << "Integral is " << Midpoint( f, a, b, N ) << '\n';
}
TheIdeasMan (6817)
@lastchance

I love when people do this: 500 LOC -> 21LOC :+D It's so good we have members with your quals and experience :+D

@Zhylian

I guess the lesson here is that not everything has to go into a class.

Don't feel bad, I think you are doing well early on in your C++ experience.
Last edited on
Zhylian (28)
Thank you guys, I feel a bit ashamed but I suppose it is part of the learning process.
TheIdeasMan (6817)
@Zhylian

Don't feel bad, everyone has to start somewhere, and everyone make mistakes.

It takes time to learn C++, there are a number of conceptual things to learn, as well as detailed things.

I can see a number good things you have done in the way you have coded your classes. I think you should give yourself a pat on the back for that :+)

C++ is a multi paradigm language:

Sometimes it is data in; process; data out. This is more like a C program but with C++ syntax. This is what lastchance has done in his example.

Sometimes using some classes is appropriate, possibly with inheritance or composition, possibly with virtual polymorphism. Design Patterns feature here, along with idioms like SOLID.

Then there is the STL. Because of the various containers, algorithms and classes in it, it is a great way of getting a lot done. A great place to start experimenting IMO.

We also have templates, which have a huge range of possibilities. Whole books written about that :+)

Finally there is Template Meta Programming (TMP) which is a whole different kettle of krawl-dads. Lots of books written about that too.

So one can tailor their code design using these paradigms - which ever is most appropriate. I guess it takes experience to realise what or which is best.



George PlusPlus (5757)
@Zhylian,

I've been self-learning C++ and the WinAPI back when Win98 was still a thing. I still make incredibly silly mistakes from time to time.
seeplus (6583)
"To err is human; to really foul things up requires a computer"

TheIdeasMan (6817)
I think what happened to the OP is one of misconception which led to a trip down a rabbit hole. The misconception happens easily, especially given other languages where everything is a class.

Sure there was a small mistake with the functions returning int rather than double, but that is a common mistake.
Zhylian (28)
Thank you guys, don't warry, I am strongly motivated to learn !
Duthomhas (13200)
Just to add some awesome history to the context here, seeplus is, in fact, properly attributing that quote by not attributing it. You can learn more here:

https://quoteinvestigator.com/2010/12/07/foul-computer/
seeplus (6583)
Updated versions could be:

"To err is human, but to really foul things up requires a robot"

or

"To err is human, but to really foul things up requires AI"
George PlusPlus (5757)
"If it's on the Internet it must be true" -- Abraham Lincoln
Topic archived. No new replies allowed.