- unsigned result = 0;
- bool errorflag = false;
- int i_num, i_den;
-
- // Check non-nested radicals (n/d)^(m/n) in ex wrapper class:
- for (int i=0; i<200; ++i) { // FIXME: run to ~200
- for (int j=2; j<13; ++j) {
- // construct an exponent 1/j...
- numeric nm(1,j);
- nm += numeric(int(20.0*rand()/(RAND_MAX+1.0))-10);
- // ...a numerator...
- do { i_num = rand(); } while (i_num == 0);
- numeric num(i_num);
- // ...and a denominator.
- do { i_den = (rand())/100; } while (i_den == 0);
- numeric den(i_den);
- // construct the radicals:
- ex radical = pow(ex(num)/ex(den),ex(nm));
- numeric floating = pow(num/den,nm);
- // test the result:
- if (is_ex_of_type(radical,numeric)) {
- clog << "(" << num << "/" << den << ")^(" << nm
- << ") should have been a product, instead it's "
- << radical << endl;
- errorflag = true;
- }
- numeric ratio = ex_to_numeric(evalf(radical))/floating;
- if (ratio>1.0001 && ratio<0.9999) {
- clog << "(" << num << "/" << den << ")^(" << nm
- << ") erroneously evaluated to " << radical;
- errorflag = true;
- }
- }
- }
- if (errorflag)
- ++result;
-
- return result;
+ unsigned result = 0;
+ bool errorflag = false;
+ int i_num, i_den;
+
+ // Check non-nested radicals (n/d)^(m/n) in ex wrapper class:
+ for (int i=0; i<200; ++i) {
+ for (int j=2; j<13; ++j) {
+ // construct an exponent 1/j...
+ numeric nm(1,j);
+ nm += numeric(int(20.0*rand()/(RAND_MAX+1.0))-10);
+ // ...a numerator...
+ do {
+ i_num = rand();
+ } while (i_num<=0);
+ numeric num(i_num);
+ // ...and a denominator.
+ do {
+ i_den = (rand())/100;
+ } while (i_den<=0);
+ numeric den(i_den);
+ // construct the radicals:
+ ex radical = pow(ex(num)/ex(den),ex(nm));
+ numeric floating = pow(num/den,nm);
+ // test the result:
+ if (is_a<numeric>(radical)) {
+ // This is very improbable with decent random numbers but it
+ // still can happen, so we better check if it is correct:
+ if (pow(radical,inverse(nm))==num/den) {
+ // Aha! We drew some lucky numbers. Nothing to see here...
+ } else {
+ clog << "(" << num << "/" << den << ")^(" << nm
+ << ") should have been a product, instead it's "
+ << radical << endl;
+ errorflag = true;
+ }
+ }
+ numeric ratio = abs(ex_to<numeric>(evalf(radical))/floating);
+ if (ratio>1.0001 && ratio<0.9999) {
+ clog << "(" << num << "/" << den << ")^(" << nm
+ << ") erroneously evaluated to " << radical;
+ errorflag = true;
+ }
+ }
+ }
+ if (errorflag)
+ ++result;
+
+ return result;