AMDG On 11/25/2011 05:11 AM, Janek Kozicki wrote:

...

...

(***) it seems strange that I need to take .value() here. The result is dimensionless, why can't it directly convert to complex<double> ? There is a risk that .value() would be expressed in millimeters or angstroms or whatever (maybe depending on dd), and then .value() would return a wrongly scaled number?

This is an idiosyncrasy of overload resolution. std::complex::operator+= is a template, so an implicit conversion doesn't work. The template argument can't be deduced.

OK. So I have to use .value().

The question here: should I expect here some scaling problems? For example:

a = 1nm b = 1km c = (a/b).value()

Will c become 1.0, with units scaling being lost, due to directly taking .value() without taking care of unit scaling. Or would it be correctly 1e-12 ?

It will be 1.0. in this case. To be safe use static_cast, or boost::implicit_cast.

In my problem I am calculating exponent then multiply and divide by length. In fact the exponent argument is supposed to be dimensionless also! And as you see I multiply wavenumber by length. I assume the conversions to ::one are good here :)

Another question:

complex<double> i(0,1); quantity<wavenumber,complex<double> > k=2*pi/lambda; // [....] a += (dd*exp(i*k*r)/r).value();

do I really need to define wavenumber k as being complex? After all it's multiplied by 'i' so it should get converted automatically?

It's the same problem as above with int/double.

...

...

quantity<length> wavelength ( 632.8 *nano*meters); quantity<length> wavelength = 632.8 *nano*meters; This isn't assignment. It's copy construction and requires an implicit conversion which is forbidden.

why it is forbidden? I think it would be more convenient to allow this.

In general, any conversion that causes the underlying value to change is explicit. In Christ, Steven Watanabe

Try D/double(steps). Matthias, I seem to recall that the reason we disabled generic multiplication was that it was causing ambiguities with user defined types. Do you think it would be a good idea to enable arithmetic with any built-in arithmetic type?

I seem to remember that we ran into some conflicts, but I can't come up with a specific example off the top of my head. I would be in favor of having arithmetic operations between quantities/units and built in types behave in the expected way.

...

(***) it seems strange that I need to take .value() here. The result is dimensionless, why can't it directly convert to complex<double> ? There is a risk that .value() would be expressed in millimeters or angstroms or whatever (maybe depending on dd), and then .value() would return a wrongly scaled number?

This is an idiosyncrasy of overload resolution. std::complex::operator+= is a template, so an implicit conversion doesn't work. The template argument can't be deduced.

Perhaps we should flesh out the toy complex implementation we wrote and incorporate it in the units library? It is a fairly common thing to use together with units and, as I understand it, the standard only defines complex for a subset of built in types...

...

Another question, why this works:

quantity<length> wavelength ( 632.8 *nano*meters);

and this does not?

quantity<length> wavelength = 632.8 *nano*meters;

Of course I know about extra copying, and waste of CPU resources on execution, when compiler isn't able to optimize this away (and usually it is able to do that). But why the heck this simple assignment can't work?

This isn't assignment. It's copy construction and requires an implicit conversion which is forbidden.

Although, in Janek's defense, this seems overly restrictive. Can we allow copy construction between scaled units and the corresponding unscaled unit without introducing undesirable side effects? Matthias

I might have got something wrong, and missed probably something in the manual... So what is a better way of writing this line in the code?

a += (dd*exp(i*k*r)/r).value();

The result of dd*exp(i*k*r)/r is dimensionless, but direct assignment (copy construction) isn't working here. And I would really admire if the code could look in a sane way, like this:

a += dd*exp(i*k*r)/r;

This would look, like I'm not really worrying about using boost::units, just benefit from all the unit checking, and no extra code to type. (and change k to use double, not complex ;) Also I am suspicious that currently my code has a bug here, exactly where the .value() is taken. Because `dd` and `r` are units of different scale.

I don't think you missed anything and I think that your second expression should work as expected. I would consider it a flaw in Boost.Units that this does not currently work and a flaw in C++ that getting a (relatively) simple mathematical concept like dimensional analysis to function in a library is as difficult as it is...

...

...

This is an idiosyncrasy of overload resolution. std::complex::operator+= is a template, so an implicit conversion doesn't work. The template argument can't be deduced.

Perhaps we should flesh out the toy complex implementation we wrote and incorporate it in the units library? It is a fairly common thing to use together with units and, as I understand it, the standard only defines complex for a subset of built in types...

that would be great. Because there is no physical reason why wavenumber k would need to be complex. It seems like allowing to perform a conversion from double to complex is reasonable. Only opposite conversion has reasons to fail. Just like converting double to int has reasons to fail. But converting int to double should work flawlessly.

I think we should be able to get such things to behave sensibly.

...

...

...

quantity<length> wavelength ( 632.8 *nano*meters); quantity<length> wavelength = 632.8 *nano*meters; Although, in Janek's defense, this seems overly restrictive. Can we allow copy construction between scaled units and the corresponding unscaled unit without introducing undesirable side effects?

I would really appreciate that :) It would simplify code. Also because I cannot put 632.8*nano*meters directly as arguments to function call.

The currently sanctioned approach to function arguments is to use templates to accept any unit of the appropriate dimension, but it does make more sense to allow implicit conversion in the special case of scaled units. The only issue is the danger of roundoff/truncation problems.

btw, how usually do you define those scaled units to be a little more convenient to use? I mean, just type: 632.8*nm Is there some header to include, to benefit from all standard SI system prefixes? It would be great, it's a standard after all :) And those short names are carefully designed to not have any unit name clashes.

I guess we could have a header corresponding to each SI unit that contains predefined scaled units for each of the SI prefixes.

Of course while writing C++ code we might end up with variables names & units names clashing, like using single letter 'm' for 'meter', but this is why we have an si:: namespace after all. And if my code grows a little bigger I will definitely switch to: 632.8*si::nm

If you will have some changes in library that will need testing I will gladly try it.

If you tell me how to define si::km, si::m, si::kN, si::GJ I can generate for you a header file with all possible units. Though 'micro' is not in latin alphabet. So micro meter is not possible, aww..

I'll try to get you something to start with later today. Matthias

Hi, Small comment on this:

'micro' is not in latin alphabet. So micro meter is not possible, aww..

You can easily use 'u' instead of mu in the greek alphabet. So micrometer would become um. That should not cause any problem right? The same would go for other quantities.. uA, uV, uN.. Best regards -ØIvind

now I started writing a simple code that would calculate a quantum wavefunction of a harmonic oscillator. It would work very quickly without units, but I want to do this using units, just to learn more using units :) And I've hit a brick wall, because I cannot get those lines to compile: quantity<reciprocal_area> lambda ( 2.0*m*E / (hbar*hbar) ); quantity<reciprocal_area> alpha ( m*omega/root<2>(hbar) ); double lambda_per_alpha ( lambda / alpha ); Even the first one has problems, to my surprise. Maybe I really should change my g++ version to a newer version? I will gladly try this if you think it makes sense. Would upgrading boost help? -- Janek Kozicki http://janek.kozicki.pl/ |

And I've hit a brick wall, because I cannot get those lines to compile:

quantity<reciprocal_area> lambda ( 2.0*m*E / (hbar*hbar) ); quantity<reciprocal_area> alpha ( m*omega/root<2>(hbar) ); double lambda_per_alpha ( lambda / alpha );

Even the first one has problems, to my surprise. Maybe I really should change my g++ 4.4.5 version to a newer version? I will gladly try this if you think it makes sense.

The first line compiles correctly using Boost 1.48; I don't know when the fix was made. The second line there is a bug/missing feature; the constants don't define pow/root operations, although I don't understand exactly why the implicit conversion to the appropriate quantity type is not working. We will fix that in the next release... For the interim, you can use root<2>(hbar.value()), although I recognize that it is somewhat confusing since, in this context, value() returns a quantity, not a value_type. However, you also have an error in the second line; alpha is not a reciprocal_area : kg * (1/s) / (kg m^2/s)^1/2 = kg^1/2 s^-1/2 m^-1. Of course, this makes the third line incorrect as well; an example of the power of Boost.Units, I guess... In general, when debugging units, it is often helpful to output to the console. Use of auto for intermediate values is also your friend... For example, here's what I did to figure out your problem: #include <iostream> #include <complex> using namespace std; //use units #include <boost/typeof/std/complex.hpp> #include <boost/units/cmath.hpp> #include <boost/units/systems/si.hpp> #include <boost/units/systems/si/prefixes.hpp> #include <boost/units/systems/si/codata/universal_constants.hpp> #include <boost/units/io.hpp> using namespace boost::units; using namespace boost::units::si; using namespace boost::units::si::constants::codata; //use units end typedef derived_dimension<length_base_dimension,-2>::type reciprocal_area_dimension; typedef unit<reciprocal_area_dimension,si::system> reciprocal_area; int main() { int n = 5; quantity<mass> m(1*kilogram); quantity<frequency> omega(1*hertz); quantity<energy> E(hbar*omega*(n+0.5)); quantity<reciprocal_area> lambda(2.0*m*E/(hbar*hbar)); std::cout << m << std::endl; std::cout << omega << std::endl; std::cout << E << std::endl; std::cout << lambda << std::endl; // std::cout << m*omega/root<2>(hbar) << std::endl; // bug - this should work std::cout << m*omega/root<2>(hbar.value()) << std::endl; // works but shouldn't be necessary and bad choice of syntax return 0; } output : 1 kg 1 s^-1 5.80014e-34 m^2 kg s^-2 1.04308e+35 m^-2 9.73782e+16 m^-1 kg^(1/2) s^(-1/2) Steven - since you implemented the constants.hpp header, I'm a little reluctant to rummage through it. Could you please consider 1) adding support for the pow<> and root<> operations? 2) deprecating constant.value_type and constant.value() in favor of constant.quantity_type and constant.quantity() or constant.mean_value() or constant.mean() or something? Also, can you remind me why we have "constant" and "physical_constant" structs? They look identical to me... Janek, I have attached headers implementing the various prefixed short cuts for the meter and joule; these two should be sufficient to extend to the rest of the SI unit system. What remains to be done is to replicate scaled_meter.hpp/scaled_length.hpp for the remaining base units in the <boost/units/base_units/si> directory and to replicate scaled_energy.hpp for all the remaining named SI units in the <boost/units/systems/si> directory. Here's my test program: #include <iostream> #include <boost/units/io.hpp> #include <boost/units/systems/si/length.hpp> #include <boost/units/systems/si/energy.hpp> #include "scaled_length.hpp" #include "scaled_energy.hpp" using namespace boost::units; using namespace boost::units::si; int main() { std::cout << 1.0*nm << "\t" << 1.5*mJ << std::endl; quantity<length> q1(1.0*nm); quantity<energy> q2(1.5*mJ); std::cout << q1 << "\t" << q2 << std::endl; auto q3(1.0*nm); auto q4(1.5*mJ); std::cout << q3 << "\t" << q4 << std::endl; return 0; } and output : 1 nm 1.5 mJ 1e-09 m 0.0015 m^2 kg s^-2 1 nm 1.5 mJ Note: if you assign these scaled quantities to unscaled quantities, the scaling information will be lost (second output line)... Matthias

Hi, thanks for examples to make scaled units (nm, km, GJ, etc..). I'll get to this soon. BTW, sorry for sending previous email twice. I pressed "send" keyboard shortcut too early. Today I wanted to finish my quantum wavefunction harmonic oscillator. I was working on it, while in parallel compiling & trying to install boost 1.48, which currently was unsuccessful (debian squeeze). BTW, if anybody has a hint about how to install boost 1.48 on debian squeeze using preferably using the "legal" `dpkg-buildpackage -rfakeroot -b` way, it would be great. But still I am quite happy today, because I managed to get a correctly working plots of this harmonic oscillator, but without units, unfortunately. One small problem is that I couldn't get to higher energies because double precision has limited factorials. I would need a working factorial of 100! which currently seems unlikely. And my currently maximum is only 25! I didn't check yet, but maybe boost has some arbitrary precision library, which would work with Boost.Units and boost::math::factorial ? So, for your pleasure I am presenting you with a harmonic oscillator code :) Unfortunately without units :) But if you will try to convert this to units, I might learn even more about using this library. Especially because dimension of wavefunction in 1 dimensional space is 1/sqrt(meter), and I have no idea how to declare this derived_dimension, this didn't work: typedef derived_dimension<length_base_dimension,root<2> >::type quantum_wavefunction_1D; In the attachment you will find a plot and code for frequency=1Hz, hbar=1 J*s, mass=1kg (plotted against a classical harmonic oscillator). Which you can compare with wikipedia and admire how similar they are :) (note that I take Planck's constant to be 1, for start). So things to do: - install boost 1.48 on my PC - make promised scaled SI units - convert this harmonic oscillator code to use Boost.Units - find some arbitrary precision library, that will work with Units and VERY big factorials, like 100! or even better 1000! I hope you enjoy this code. I like it :) best regards -- Janek Kozicki http://janek.kozicki.pl/ |

Paul A. Bristow said: (by the date of Wed, 30 Nov 2011 11:08:53 -0000)

...

- find some arbitrary precision library, that will work with Units and VERY big factorials, like 100! or even better 1000!

There are indeed at least two related efforts on arbitrary precision library going on.

A multiprecision library by Christopher Kormanyos based on his e_float library which is in the sandbox .

and, based on this, an extension to Boost.Math (already user-defined type ready and able to use NTL and GMP) that John Maddock is currently busy working on.

These should allow you to calculate factorial 100! to a few hundred decimal digits.

For your amusement, 100! is (to 50 decimal digits)

100! = 9.332621544394415268169923885626670049071596826438162146859e+157

I'm using apcalc, a very nice software ;-) 933262154439441526816992388562667004907159682643816214685929638952175999932\ 29915608941463976156518286253697920827223758251185210916864000000000000000000000000

But I'm not at all sure that this is the way forward for your purpose. For example, in calculating gamma functions, it is customary to use logs to avoid overflow.

And we need to make sure we can walk before we can run - combining arbitrary precision with Boost.Units (where double precision is usually enough to represent the accuracy of physical things).

hmm, you are right. I probably will need to go into logs. Boost.Units has no effect on runtime, it checks units during compilation. So using them shouldn't affect at all. I'm glad that I installed 1.48 with Steven's help. And now I can convert my code to use units :) first thing will be to figure out how to declare a type meter^-0.5 maybe my previous try that didn't work with boost 1.42 will start to work with 1.48: typedef derived_dimension<length_base_dimension,root<2> >::type quantum_wavefunction_1D; best regards -- Janek Kozicki http://janek.kozicki.pl/

Matthias Schabel said: (by the date of Wed, 30 Nov 2011 07:51:19 -0800)

...

first thing will be to figure out how to declare a type meter^-0.5 maybe my previous try that didn't work with boost 1.42 will start to work with 1.48:

typedef derived_dimension<length_base_dimension,root<2> >::type quantum_wavefunction_1D;

I tried to send a working version of your code yesterday, but the message was bounced because the attachments were too big; I will resend when I get to the office. For now, your problem here is twofold:

1) derived_dimension is a convenience typedef for integer powers; you need to use make_dimension_list<> 2) once you have defined the dimension, you also need to typedef a unit in the SI system corresponding to your desired length^(-1/2) dimension:

typedef unit<quantum_wavefunction_1D,si::system> quantum_wavefunction_1D_unit

great! Thanks for hint and your effort to make it work. I'm not sure yet how to use make_dimension_list<> so perhaps better if I wait for your code, than spend time to figure it out :) Fact that it is one-dimensional wavefunction is important. Because in 3 dimensional space quantum wavefunction has dimension of (length)^(1.5) == square root of volume. The reasoning is following. Let's say that psi(r) is a wavefunction (where 'r' are space coordinates). Then it's squared magnitude is probability. And total probability in whole infinite space must be 1.0 dimensionless, so: integral abs(psi(r))^2 dr = 1.0 dimensionless this is why quantum wavefunction depending on number of dimensions is 1D: (length)^(-0.5) ; 2D: (length)^(-1.0) ; 3D: (length)^(-1.5) It would be great if you added this to the directory `physical_dimensions`. Also for 2D and 3D. Now I am not sure about naming. I used "1D" name here. If you have some other convention instead of "1D" already in the naming, then use that convention. Another thing to note: this is quantum wavefunction in spatial representation, it's on *length*. It's possible to have other representations as well, for example on *momentum*. Then in 1D it would be (momentum)^0.5 http://en.wikipedia.org/wiki/Momentum_space There are a dozen of other representations possible. They are not very popular though. So I don't know, maybe it would be best if you added 6 wavefunction types: 1D, 2D, 3D for length and for momentum. I'm not sure about their naming though. If one says 'wavefunction' it is assumed that it is in spatial coordinates. If it's in momentum space, then it's said explicitly 'wavefunction in momentum space' best regards -- Janek Kozicki http://janek.kozicki.pl/ |

Here is a corrected version of your code with a few notes: #include <iostream> #include <complex> #include <boost/foreach.hpp> using namespace std; // use rational numbers #include <boost/rational.hpp> // use math constants #include <boost/math/constants/constants.hpp> // use math factorial #include <boost/math/special_functions/factorials.hpp> //using namespace boost::math // use units #include <boost/typeof/std/complex.hpp> // MCS - need to include this header for make_dimension_list #include <boost/mpl/list.hpp> #include <boost/units/cmath.hpp> // MCS - need to include this header #include <boost/units/pow.hpp> #include <boost/units/io.hpp> #include <boost/units/systems/si.hpp> #include <boost/units/systems/si/prefixes.hpp> #include <boost/units/systems/si/codata/electron_constants.hpp> #include <boost/units/systems/si/codata/physico-chemical_constants.hpp> #include <boost/units/systems/si/codata/universal_constants.hpp> using namespace boost::units; using namespace boost::units::si; // MCS - need codata namespace to find hbar using namespace boost::units::si::constants::codata; //use units end vector<boost::rational<signed long> > hermite_polynomial_coefficients(int order, boost::rational<signed long> lambda_per_alpha) { vector<boost::rational<signed long> > c_even; // can't use static, because lambda_per_alpha changes between calls vector<boost::rational<signed long> > c_odd; if(c_even.size() == 0) { c_even.push_back(1); c_even.push_back(0); } if(c_odd.size() == 0) { c_odd.push_back(0); c_odd.push_back(1); } vector<boost::rational<signed long> >& a( (order % 2 == 0) ?c_even:c_odd ); if(order+1 < a.size()) return a; for(int i = a.size() ; i <= order ; ++i) { boost::rational<signed long> next = a[i-2]*(2*(i-2)+1-lambda_per_alpha)/(i*(i-1)); a.push_back(next); } // BOOST_FOREACH(boost::rational<signed long> v,a) // std::cout << v << " "; // std::cout << "\n"; return a; } vector<boost::rational<signed long> > hermite_polynomial_scaled(int order, boost::rational<signed long> lambda_per_alpha) { vector<boost::rational<signed long> > a(hermite_polynomial_coefficients(order,lambda_per_alpha)); vector<boost::rational<signed long> > b; boost::rational<signed long> factor(std::pow(2.,order)/a[order]); for(int i = 0 ; i<=order ; ++i) b.push_back(a[i]*factor); return b; } // MCS - derived_dimension only works with integer powers; use make_dimension_list here instead typedef make_dimension_list<boost::mpl::list<dim<length_base_dimension,static_rational<-1,2> > > >::type quantum_wavefunction_1D_dimension; typedef unit<quantum_wavefunction_1D_dimension,si::system> quantum_wavefunction_1D; // this function returns 1/sqrt(meter) quantity<quantum_wavefunction_1D,complex<double> > quantum_oscillator_wavefunction( int n // n - order of wavefunction , const quantity<length> x // position , const quantity<mass> mass // mass - mass of particle , const quantity<frequency> omega // oscillator frequency ) { // MCS - need to define unit typedef in addition to dimension typedef derived_dimension<length_base_dimension,-2>::type reciprocal_area_dimension; typedef unit<reciprocal_area_dimension,si::system> reciprocal_area; quantity<energy_time> h_(hbar); quantity<energy> E(h_*omega*(n+0.5)); // energy of oscillator quantity<reciprocal_area> lambda(2.0*mass*E/(h_*h_)), alpha(mass*omega/h_); quantity<wavenumber> alpha_root(root<2>(alpha)); double lambda_per_alpha(2*n+1); const double pi(boost::math::constants::pi<double>()), fac(boost::math::factorial<double>(n)); // N_n should be 1/sqrt(meter) because it has a dimension of quantum_wavefunction_1D const quantity<quantum_wavefunction_1D,complex<double> > N_n(root<2>(root<2>(alpha/pi)/(std::pow(2.,n)*fac))); const complex<double> i(0,1); quantity<quantum_wavefunction_1D,complex<double> > result(0/root<2>(meter)); vector<boost::rational<signed long> > a(hermite_polynomial_scaled(n,lambda_per_alpha)); for(int v=0 ; v<=n ; v++) { // result += N_n*a[v]*pow<v>(alpha_root*x)*exp(-0.5*alpha*pow<2>(x)); // MCS - a[v] is boost::rational -> need to convert to double via rational_cast // MCS - pow<v> does not work because v must be known at compile time const double a_v(boost::rational_cast<double>(a[v])), tmp1(a_v*pow(alpha_root*x,v)); const complex<double> tmp2(exp(-0.5*alpha*pow<2>(x))), tmp3(tmp1*tmp2); result += tmp3*N_n; } return result; } int main() { // MCS - change -1 -> +1 here otherwise loop doesn't run quantity<length> s(15*nano*meter), ds(0.01*nano*meter); std::cout << s << "\t" << ds << std::endl; for (quantity<length> x=-s ; x<=s ; x+=ds) { const quantity<temperature> T(273.15*kelvin); const quantity<frequency> omega(k_B*T/hbar.value()); std::cout << x << " " << quantum_oscillator_wavefunction(25,x,m_e,omega) << "\n"; } for (quantity<length> x=-s ; x<=s ; x+=ds) { const quantity<temperature> T(273.15*kelvin); const quantity<frequency> omega(k_B*T/hbar.value()); const quantity<quantum_wavefunction_1D,complex<double> > psi(quantum_oscillator_wavefunction(25,x,m_e,omega)); std::cout << x.value() << " " << psi.value().real() << "\t" << psi.value().imag() << "\n"; } return 0; } and the n=25 wavefunction: Cheers, Matthias