开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 查看源码
_real.c
资源名称：leda.tar.gz　[点击查看]
上传用户：gzelex
上传日期：2007-01-07
资源大小：707k
文件大小：16k
源码类别：

数值算法/人工智能

开发平台：

MultiPlatform

_real.c：源码内容
  1. /*******************************************************************************
  2. +
  3. +  LEDA-R  3.2.3
  4. +
  5. +  _real.c
  6. +
  7. +  Copyright (c) 1995  by  Max-Planck-Institut fuer Informatik
  8. +  Im Stadtwald, 66123 Saarbruecken, Germany     
  9. +  All rights reserved.
  10. + 
  11. *******************************************************************************/
  14. /* This file has been automatically generated from "real.w"
  15.    by CTANGLE (Version 3.1 [km2c]) */
  16. /*12:*/
  18. #include <LEDA/real.h>
  19. /*16:*/
  21. #include <math.h>
  25. /*:16*/
  27. /*14:*/
  29. enum constructortype{DOUBLE= 0,BIGFLOAT= 1,NEGATION= 2,
  30. SQUAREROOT= 3,ADDITION= 4,SUBTRACTION= 5,MULTIPLICATION= 6,
  31. DIVISION= 7};
  33. struct app_bf
  34. {bigfloat NUM_BF;
  35. bigfloat ERROR_BF;
  36. app_bf(){}
  37. app_bf(const bigfloat&x,const bigfloat&error_x)
  38. {NUM_BF= x;ERROR_BF= error_x;}
  39. };
  41. /*:14*/
  43. /*3:*/
  45. const double eps= ldexp(1,-53);/* $eps = 2^{-53}$ */
  46. const double correction= 1+8*eps;
  47. const double MaxError= ldexp(1,900);/* $2^{900}$ */
  48. const double MinDbl= ldexp(1,-1022);
  49. const double MaxDbl= ldexp(2-ldexp(1,-52),1023);
  51. /*:3*//*4:*/
  53. const bigfloat zero= double(0);
  54. const bigfloat one= double(0);
  55. const bigfloat two= double(2);
  56. const bigfloat onehalf= double(1)/double(2);
  57. const bigfloat onefourth= double(1)/double(4);
  58. const bigfloat oneeighth= double(1)/double(8);
  59. const bigfloat eps_prec= eps;
  60. const bigfloat correction_bf= correction;
  62. /*:4*/
  64. /*15:*/
  66. class real_rep
  67. {
  68. friend class real;
  69. /*9:*/
  71. friend real operator+(const real&,const real&);
  72. /*{Mbinop addition
  73. }*/
  74. friend real operator-(const real&,const real&);
  75. /*{Mbinop subtraction
  76. }*/
  77. friend real operator*(const real&,const real&);
  78. /*{Mbinop multiplication
  79. }*/
  80. friend real operator/(const real&,const real&);
  81. /*{Mbinop division
  82. }*/
  83. friend real operator-(const real&);
  84. /*{Munop negation
  85. }*/
  86. friend real sqrt(const real&);
  87. /*{Mop squareroot operation
  88. }*/
  90. /*:9*/
  93. double NUM;
  94. double ERROR;
  95. app_bf*APP_BF_PTR;
  96. constructortype CON;
  97. real_rep*OP1,*OP2;
  98. int count;
  100. real_rep(double= 0);
  101. real_rep(const bigfloat&);
  102. real_rep(const bigrational&);
  104. ~real_rep();
  106. LEDA_MEMORY(real_rep)
  108. inline bigfloat&num_bf() const {return APP_BF_PTR->NUM_BF;}
  109. inline bigfloat&error_bf() const {return APP_BF_PTR->ERROR_BF;}
  110. void init_app_bf();
  111. void adjust_dbl();
  113. int sign();
  114. int sign(const bigfloat&);
  115. inline void improve(const bigfloat&);
  116. inline void improve(long k);
  117. inline void compute(long);
  119. /*43:*/
  121. friend void Mignotte_parameters(const real_rep&,integer&,integer&);
  123. /*:43*/
  125. };
  127. /*:15*/
  129. // INSERTED 
  130. GenPtr real::copy() const { PTR->count++; return PTR; }
  131. void   real::clear()      { if (PTR && --(PTR->count)==0) delete PTR; }
  134. /*17:*/
  136. real::real(){PTR= new real_rep(0);}
  137. real::real(double x){PTR= new real_rep(x);}
  138. real::real(int x){PTR= new real_rep(double(x));}
  139. real::real(const bigrational&x){PTR= new real_rep(x);}
  140. real::real(const bigfloat&x){PTR= new real_rep(x);}
  141. real::real(const integer&x){PTR= new real_rep(bigfloat(x));}
  142. real::real(real_rep&x_rep){PTR= &x_rep;}
  143. real::real(const real&x){PTR= x.PTR;x.PTR->count++;}
  144. real&real::operator= (const real&x)
  145. {x.PTR->count++;
  146. if(PTR && --PTR->count==0)delete PTR;
  147. PTR= x.PTR;
  148. return(*this);
  149. }
  151. /*:17*//*18:*/
  153. real_rep::real_rep(double x)
  154. {APP_BF_PTR= nil;count= 1;
  155. OP1= OP2= nil;
  156. NUM= x;ERROR= 0;CON= DOUBLE;
  157. }
  158. real_rep::real_rep(const bigfloat&x)
  159. {APP_BF_PTR= new app_bf(x,0);count= 1;
  160. OP1= nil;OP2= nil;CON= BIGFLOAT;
  161. adjust_dbl();
  162. }
  163. real_rep::real_rep(const bigrational&x)
  164. {
  165. APP_BF_PTR= nil;count= 1;
  166. OP1= new real_rep(bigfloat(x.numerator()));
  167. OP2= new real_rep(bigfloat(x.denominator()));
  168. CON= DIVISION;
  169. adjust_dbl();
  170. }
  172. /*:18*//*19:*/
  174. real::~real()
  175. {PTR->count--;
  176. if(PTR->count==0)
  177. delete PTR;
  178. }
  180. real_rep::~real_rep()
  181. {
  182. if(APP_BF_PTR!=nil)
  183. delete APP_BF_PTR;
  184. if((OP1!=nil)&&(--OP1->count==0))
  185. delete OP1;
  186. if((OP2!=nil)&&(--OP2->count==0))
  187. delete OP2;
  188. #ifdef DEBUG_REAL
  189. APP_BF_PTR= nil;
  190. OP1= OP2= nil;
  191. #endif
  192. }
  195. /*:19*/
  197. /*20:*/
  199. void real_rep::init_app_bf()
  200. {
  201. if(APP_BF_PTR==nil)
  202. if(ERROR!=Infinity)
  203. APP_BF_PTR= 
  204. new app_bf(bigfloat(NUM),fabs(bigfloat(ERROR*NUM*correction)));
  205. else
  206. compute(52);
  207. }
  210. /*:20*//*21:*/
  212. void real_rep::adjust_dbl()
  213. {
  214. if((fabs(num_bf())>bigfloat(MaxDbl))
  215. ||(error_bf()>bigfloat(MaxDbl))
  216. )
  217. {
  218. ERROR= Infinity;NUM= 1;return;
  219. }
  220. if((num_bf()==zero)&&(error_bf()==zero))
  221. {
  222. NUM= ERROR= 0;return;
  223. }
  224. double n_head= num_bf().todouble();
  225. double e_head= error_bf().todouble();
  226. if(n_head!=0)
  227. {NUM= n_head;
  228. ERROR= eps+Max(MinDbl,e_head)/fabs(n_head);
  229. }
  230. else
  231. {NUM= e_head+2*MinDbl;
  232. ERROR= 2;
  233. }
  234. ERROR= ERROR*correction;
  235. if(ERROR>=MaxError)
  236. {NUM= 1;ERROR= Infinity;
  237. }
  238. }
  240. /*:21*//*22:*/
  243. double real::todouble()const
  244. {return PTR->NUM;}
  245. double real::get_double_error()const
  246. {return PTR->ERROR;}
  247. bigfloat real::tobigfloat()const
  248. {PTR->init_app_bf();return PTR->num_bf();}
  249. integer real::get_precision()const
  250. {PTR->init_app_bf();return-PTR->error_bf().get_exponent();}
  252. ostream&operator<<(ostream&out, const real&x)
  253. {out<<x.tobigfloat();return out;}
  255. istream&operator>>(istream&in,real&x)
  256. {
  257. double x_num;in>>x_num;
  258. x= real(x_num);return in;
  259. }
  263. /*:22*/
  265. /*23:*/
  267. real operator+(const real&x,const real&y)
  268. {
  269. real_rep&z_rep= *new real_rep();
  270. /* this represents the new |real| element |z| that we return later */
  272. z_rep.CON= ADDITION;
  273. z_rep.OP1= x.PTR;z_rep.OP2= y.PTR;
  274. z_rep.OP1->count++;z_rep.OP2->count++;
  275. /* the new element |z| now points to |x| and |y|, so their counters have
  276.        to be increased by one
  277.     */
  279. double&x_num= x.PTR->NUM;
  280. double&y_num= y.PTR->NUM;
  281. double&NUM= z_rep.NUM;
  282. NUM= x_num+y_num;/* approximation of $x+y$ */
  284. if(fabs(NUM)>MinDbl)
  285. {/* regular case */
  286. z_rep.ERROR= fabs(x_num/NUM)*x.PTR->ERROR
  287. +fabs(y_num/NUM)*y.PTR->ERROR
  288. +eps;
  289. }
  290. else
  291. {/* underflow is possible */
  293. NUM= fabs(x_num)*x.PTR->ERROR+fabs(y_num)*y.PTR->ERROR+4*MinDbl;
  294. z_rep.ERROR= 2;
  295. }
  297. z_rep.ERROR*= correction;
  299. if(!((fabs(NUM)+x.PTR->ERROR+y.PTR->ERROR<Infinity)
  300. &&(z_rep.ERROR<MaxError)
  301. )
  302. )
  303. /* if one of the quantities fabs(NUM), x.PTR->ERROR, y.PTR->ERROR, and
  304.         z_rep.ERROR is either Infinity or NaN, we set the error of the
  305.         result to Infinity
  306.      */
  307. z_rep.ERROR= Infinity;
  309. return real(z_rep);
  310. }
  312. real operator-(const real&x,const real&y)
  313. {
  314. real_rep&z_rep= *new real_rep();
  315. z_rep.CON= SUBTRACTION;
  316. z_rep.OP1= x.PTR;z_rep.OP2= y.PTR;
  317. z_rep.OP1->count++;z_rep.OP2->count++;
  318. if((x.PTR->ERROR==Infinity)||(y.PTR->ERROR==Infinity))
  319. {z_rep.ERROR= Infinity;return real(z_rep);
  320. }
  321. double x_num= x.PTR->NUM,y_num= y.PTR->NUM;
  322. double&NUM= z_rep.NUM;
  323. NUM= x_num-y_num;/* approximation of $x-y$ */
  324. if(fabs(NUM)<=MinDbl)
  325. {
  326. NUM= fabs(x_num)*x.PTR->ERROR+fabs(y_num)*y.PTR->ERROR+4*MinDbl;
  327. z_rep.ERROR= 2;
  328. }
  329. else
  330. {z_rep.ERROR= fabs(x_num/NUM)*x.PTR->ERROR
  331. +fabs(y_num/NUM)*y.PTR->ERROR
  332. +eps;
  333. }
  334. z_rep.ERROR*= correction;
  336. if(!((fabs(NUM)+x.PTR->ERROR+y.PTR->ERROR<Infinity)
  337. &&(z_rep.ERROR<MaxError)
  338. )
  339. )
  340. z_rep.ERROR= Infinity;
  342. return real(z_rep);
  343. }
  345. /*:23*//*24:*/
  347. real operator*(const real&x,const real&y)
  348. {
  349. real_rep&z_rep= *new real_rep();
  350. z_rep.CON= MULTIPLICATION;
  351. z_rep.OP1= x.PTR;z_rep.OP2= y.PTR;
  352. z_rep.OP1->count++;z_rep.OP2->count++;
  353. double&x_num= x.PTR->NUM;
  354. double&y_num= y.PTR->NUM;
  355. double&qx= x.PTR->ERROR;
  356. double&qy= y.PTR->ERROR;
  357. double&NUM= z_rep.NUM;
  358. NUM= x_num*y_num;
  360. if(!(fabs(NUM)<MinDbl))
  361. {
  362. z_rep.ERROR= qx+qy+qx*qy+eps;
  363. }
  364. else
  365. {if((x_num!=0)&&(y_num!=0))
  366. {/* underflow occurred in the computation of |z_rep.NUM| */
  367. z_rep.NUM= MinDbl*(1+qx)*(1+qy);
  368. z_rep.ERROR= 2;
  369. }
  370. else
  371. z_rep.ERROR= 0;
  372. }
  373. z_rep.ERROR*= correction;
  375. if(!((fabs(NUM)+qx+qy<Infinity)&&(z_rep.ERROR<MaxError)))
  376. z_rep.ERROR= Infinity;
  378. return real(z_rep);
  379. }
  381. /*:24*//*25:*/
  383. real operator/(const real&x,const real&y)
  384. {
  385. real_rep&z_rep= *new real_rep();
  386. z_rep.CON= DIVISION;
  387. z_rep.OP1= x.PTR;z_rep.OP2= y.PTR;
  388. z_rep.OP1->count++;z_rep.OP2->count++;
  389. double x_num= x.PTR->NUM,y_num= y.PTR->NUM;
  390. double&NUM= z_rep.NUM;
  391. if(!(y.PTR->ERROR<1))
  392. {z_rep.ERROR= Infinity;return real(z_rep);
  393. }
  394. if(y_num==0)
  395. error_handler(1,"real::operator/:Division by zero");
  396. NUM= x_num/y_num;
  397. if(fabs(NUM)>=MinDbl)
  398. z_rep.ERROR= (x.PTR->ERROR+y.PTR->ERROR)/(1-y.PTR->ERROR)
  399. +eps;
  400. else
  401. {
  402. if(x_num!=0)
  403. {
  404. NUM= MinDbl*(1+x.PTR->ERROR)/(1-y.PTR->ERROR);
  405. z_rep.ERROR= 2;
  406. }
  407. else z_rep.ERROR= 0;
  408. }
  409. z_rep.ERROR*= correction;
  411. if(!((fabs(NUM)+x.PTR->ERROR!=Infinity)||(z_rep.ERROR<MaxError)))
  412. z_rep.ERROR= Infinity;
  413. return real(z_rep);
  414. }
  416. /*:25*//*26:*/
  418. real sqrt(const real&x)
  419. {real_rep&z_rep= *new real_rep();
  420. z_rep.CON= SQUAREROOT;
  421. z_rep.OP1= x.PTR;z_rep.OP1->count++;
  422. z_rep.OP2= nil;
  423. double x_num= x.PTR->NUM;
  424. if(x.PTR->ERROR>=1)
  425. {z_rep.ERROR= Infinity;return real(z_rep);
  426. }
  427. if(x_num==0)
  428. {z_rep.NUM= z_rep.ERROR= 0;return real(z_rep);
  429. }
  430. if(x_num<0)
  431. error_handler(-1,"real::sqrt: sqrt of negative number");
  432. z_rep.NUM= sqrt(x_num);
  433. z_rep.ERROR= x.PTR->ERROR+eps;
  434. z_rep.ERROR*= correction;
  435. return real(z_rep);
  436. }
  438. /*:26*//*27:*/
  440. real operator-(const real&x)
  441. {real_rep&z_rep= *new real_rep();
  442. z_rep.CON= NEGATION;
  443. z_rep.OP1= x.PTR;z_rep.OP1->count++;z_rep.OP2= nil;
  444. z_rep.NUM= -x.PTR->NUM;z_rep.ERROR= x.PTR->ERROR;
  445. return real(z_rep);
  446. }
  448. /*:27*//*28:*/
  450. real operator+= (real&x,real&y){x= x+y;return x;}
  451. real operator-= (real&x,real&y){x= x-y;return x;}
  452. real operator*= (real&x,real&y){x= x*y;return x;}
  456. /*:28*/
  458. /*32:*/
  460. long needed_places(const bigfloat&epsilon)
  461. {return-epsilon.get_exponent().tolong();}
  463. /*:32*/
  465. /*29:*/
  467. void real::improve(const integer&prec){PTR->improve(pow2(-prec));}
  468. void real::improve(long p){PTR->improve(pow2(-p));}
  469. void real::compute_in(long p){PTR->compute(p);}
  470. void real::compute_up_to(long p)
  471. {PTR->init_app_bf();
  472. if(PTR->sign()==0)return;
  473. bigfloat q= (fabs(tobigfloat())+pow2(-get_precision()))*bigfloat(pow2(-p-1));
  474. PTR->improve(q);
  475. PTR->APP_BF_PTR->NUM_BF.round(p);
  476. PTR->APP_BF_PTR->ERROR_BF= q;
  477. PTR->adjust_dbl();
  478. }
  480. /*:29*//*30:*/
  482. void real_rep::improve(const bigfloat&q)
  483. {
  484. if(q<=zero)
  485. error_handler(-1,"real_rep::improve: quality <= 0");
  486. int signum;
  487. long p;
  488. bigfloat needed_precision,
  489. easy_value,
  490. qx,qy,
  491. y_low,
  492. x_high,x_low,
  493. qsquare;
  494. init_app_bf();
  495. if(error_bf()<=q)return;
  497. switch(CON)
  498. {
  499. case DOUBLE:case BIGFLOAT:return;
  500. case NEGATION:
  501. OP1->improve(q);
  502. num_bf()= -OP1->num_bf();
  503. break;
  504. case ADDITION:/*33:*/
  506. OP1->improve(onefourth*q);
  507. OP2->improve(onefourth*q);
  508. easy_value= OP1->num_bf()+OP2->num_bf();
  509. if(easy_value==zero)
  510. num_bf()= zero;
  511. else
  512. {
  513. needed_precision= onefourth*q/fabs(easy_value);
  514. p= needed_places(needed_precision);
  515. num_bf()= add(OP1->num_bf(),OP2->num_bf(),p);
  516. }
  517. /*:33*/
  518. ;break;
  519. case SUBTRACTION:/*34:*/
  521. OP1->improve(onefourth*q);
  522. OP2->improve(onefourth*q);
  523. easy_value= OP1->num_bf()+OP2->num_bf();
  524. if(easy_value==zero)
  525. num_bf()= 0;
  526. else
  527. {
  528. needed_precision= onefourth*q/fabs(easy_value);
  529. p= needed_places(needed_precision);
  530. num_bf()= sub(OP1->num_bf(),OP2->num_bf(),p);
  531. }
  533. /*:34*/
  534. ;break;
  535. case MULTIPLICATION:/*35:*/
  537. OP2->init_app_bf();
  538. qx= q*oneeighth/(fabs(OP2->num_bf())+OP2->error_bf());
  539. OP1->improve(qx);
  540. qy= q*oneeighth/fabs(OP1->num_bf());
  541. OP2->improve(qy);
  542. easy_value= OP1->num_bf()*OP2->num_bf();
  543. if(easy_value==zero)
  544. num_bf()= zero;
  545. else
  546. {
  547. needed_precision= q*onefourth/fabs(easy_value);
  548. p= needed_places(needed_precision);
  549. num_bf()= mul(OP1->num_bf(),OP2->num_bf(),p);
  550. }
  553. /*:35*/
  554. ;break;
  555. case DIVISION:/*36:*/
  557. signum= OP2->sign();
  558. if(signum==0)error_handler(0,"real::improve: division by zero");
  559. OP2->init_app_bf();
  560. y_low= fabs(OP2->num_bf())-error_bf();
  561. qx= q*y_low*oneeighth;
  562. OP1->improve(qx);
  563. if(OP1->num_bf()==zero)
  564. num_bf()= 0;
  565. else
  566. {
  567. qy= q*y_low*y_low*oneeighth/fabs(OP1->num_bf());
  568. OP2->improve(qy);
  569. needed_precision= q*onefourth*fabs(OP2->num_bf()/OP1->num_bf());
  570. p= needed_places(needed_precision);
  571. num_bf()= div(OP1->num_bf(),OP2->num_bf(),p);
  572. }
  574. /*:36*/
  575. ;break;
  576. case SQUAREROOT:/*37:*/
  578. OP1->init_app_bf();
  579. x_low= OP1->num_bf()-OP1->error_bf();
  580. x_high= OP1->num_bf()+OP1->error_bf();
  581. qsquare= q*q;
  583. if((x_low<=onefourth*qsquare)&&(x_high>=onehalf*qsquare))
  584. {
  585. OP1->improve(oneeighth*qsquare);
  586. x_low= OP1->num_bf()-OP1->error_bf();
  587. x_high= OP1->num_bf()+OP1->error_bf();
  588. }
  589. if(x_high<zero)
  590. error_handler(-1,"real::improve: sqrt needs an argument >= 0");
  591. if(x_low>onefourth*qsquare)
  592. {
  593. qx= onefourth*q*sqrt(x_low);
  594. OP1->improve(qx);
  595. needed_precision= onefourth*q/sqrt(x_high);
  596. p= needed_places(needed_precision);
  597. num_bf()= sqrt(fabs(OP1->num_bf()),p);
  598. }
  599. else
  600. num_bf()= 0;
  603. /*:37*/
  604. ;break;
  605. }
  606. error_bf() = q;
  607. adjust_dbl();
  608. }
  610. /*:30*//*38:*/
  612. void real_rep::compute(long p)
  613. {
  614. if(p<=0)
  615. error_handler(1,"real_rep::compute needs an argument > 0");
  616. bigfloat*x_ptr= nil,*qx_ptr= nil,*y_ptr= nil,*qy_ptr= nil;
  617. if(OP1!=nil)
  618. {OP1->compute(p);
  619. x_ptr= &(OP1->num_bf());
  620. qx_ptr= &(OP1->error_bf());
  621. }
  622. if(OP2!=nil)
  623. {OP2->compute(p);
  624. y_ptr= &(OP2->num_bf());
  625. qy_ptr= &(OP2->error_bf());
  626. }
  627. bigfloat&x= *x_ptr,&y= *y_ptr,&qx= *qx_ptr,&qy= *qy_ptr,
  628. x_low,x_high,y_low;
  629. if(APP_BF_PTR==nil)APP_BF_PTR= new app_bf();
  631. int signum;
  632. bigfloat prec= pow2(-p);
  634. switch(CON)
  635. {
  636. case DOUBLE:
  637. num_bf()= bigfloat(NUM);
  638. error_bf()= 0;break;
  639. case BIGFLOAT:break;
  640. case NEGATION:
  641. num_bf()= -x;
  642. error_bf()= qx;break;
  643. case ADDITION:
  644. num_bf()= add(x,y,p);
  645. error_bf()= fabs(num_bf())*prec+qx+qy;break;
  646. case SUBTRACTION:
  647. num_bf()= sub(x,y,p);
  648. error_bf()= fabs(num_bf())*prec+qx+qy;break;
  649. case MULTIPLICATION:
  650. num_bf()= mul(x,y,p);
  651. error_bf()= fabs(num_bf())*prec+fabs(x)*qy+fabs(y)*qx+qx*qy;
  652. break;
  653. case DIVISION:
  654. if(qy>=fabs(y))
  655. {signum= OP2->sign();
  656. if(signum==0)error_handler(-2,"real::compute: division by zero");
  657. }
  658. y_low= fabs(OP2->num_bf())-OP2->error_bf();
  659. num_bf()= div(x,y,p);
  660. error_bf()= fabs(num_bf())*prec+(qy*fabs(num_bf())+qx)/y_low;
  661. break;
  662. case SQUAREROOT:
  663. x_low= x-qx;x_high= x+qx;
  664. if(x_high<zero)
  665. error_handler(-1,"real::compute: squareroot needs argument >= 0");
  666. if(x_low>zero)
  667. {
  668. num_bf()= sqrt(x,p);
  669. error_bf()= num_bf()*prec+qx/sqrt(x_low);
  670. }
  671. else
  672. {
  673. num_bf()= sqrt(qx,p);
  674. error_bf()= two*num_bf();
  675. }
  676. }
  677. error_bf()= error_bf()*bigfloat(correction);
  678. adjust_dbl();
  679. }
  682. /*:38*/
  684. /*39:*/
  686. int sign(double x)
  687. {if(x==0)return 0;else return(x>0?1:-1);}
  689. /*:39*//*40:*/
  691. int real::sign() {return PTR->sign(0);}
  692. int real::sign(const integer&prec) {return PTR->sign(pow2(-prec));}
  693. int real::sign(long p) {return PTR->sign(pow2(-p));}
  694. int real_rep::sign() {return sign(0);}
  695. int real_rep::sign(const bigfloat&q_initial)
  696. {integer M= 0,deg= 1;
  697. bigfloat q= q_initial;
  698. /*41:*/
  700. if((ERROR!=Infinity)&&((NUM==0)||(ERROR<1)))
  701. return::sign(NUM);
  702. if((APP_BF_PTR!=nil)&&
  703. ((error_bf()==zero)||(error_bf()<fabs(num_bf())))
  704. )
  705. return::sign(num_bf());
  707. /*:41*/
  709. if(q==zero)/*45:*/
  711. Mignotte_parameters(*this,M,deg);
  713. q= pow2(-M)*correction_bf;
  715. /*:45*/
  717. /*42:*/
  719. bigfloat current_precision= eps_prec;
  720. init_app_bf();
  721. if(error_bf()<one)
  722. current_precision= error_bf()*eps_prec;
  723. improve(current_precision);
  724. while((current_precision>q)&&(error_bf()>=fabs(num_bf())))
  725. {current_precision= current_precision*current_precision;
  726. improve(current_precision);
  727. }
  728. if(fabs(num_bf())>error_bf())
  729. return::sign(num_bf());
  730. return 0;
  732. /*:42*/
  734. }
  736. /*:40*//*44:*/
  738. void Mignotte_parameters(const real_rep&x,integer&M,integer&deg)
  739. {
  740. long exp,M_int;
  741. integer M1,M2,deg1,deg2;
  742. switch(x.CON)
  743. {case DOUBLE:
  744. deg= 1;
  745. exp= bigfloat(x.NUM).get_exponent().tolong();
  746. if(x.NUM==int(x.NUM))/* x.NUM is an integer */
  747. {M= exp+1;break;}
  748. if(exp>=52)
  749. {M_int= exp+1;M= M_int;break;}
  750. if(exp<0){M_int= 52-exp;M= M_int;break;}
  751. M= 53;break;
  752. case BIGFLOAT:if(x.APP_BF_PTR->NUM_BF!=zero)
  753. M= (x.APP_BF_PTR->NUM_BF).get_exponent();
  754. else
  755. M= 0;
  756. deg= 1;break;
  757. case NEGATION:Mignotte_parameters(*x.OP1,M,deg);break;
  758. case MULTIPLICATION:case DIVISION:
  759. Mignotte_parameters(*x.OP1,M1,deg1);
  760. Mignotte_parameters(*x.OP2,M2,deg2);
  761. M= M1*deg2+M2*deg1;
  762. deg= deg1*deg2;
  763. break;
  764. case ADDITION:case SUBTRACTION:
  765. Mignotte_parameters(*x.OP1,M1,deg1);
  766. Mignotte_parameters(*x.OP2,M2,deg2);
  767. M= deg1*deg2+M1*deg2+M2*deg1;
  768. deg= deg1*deg2;
  769. break;
  770. case SQUAREROOT:
  771. Mignotte_parameters(*x.OP1,M1,deg1);
  772. M= M1;deg= deg1<<1;
  773. break;
  774. }
  775. }
  777. /*:44*//*46:*/
  779. int operator==(const real&x,const real&y)
  780. {return((y-x).sign()==0?1:0);}
  781. int operator<(const real&x,const real&y)
  782. {return((y-x).sign()>0?1:0);}
  783. int operator>(const real&x,const real&y)
  784. {return((y-x).sign()<0?1:0);}
  785. int operator<=(const real&x,const real&y)
  786. {return((y-x).sign()>=0?1:0);}
  787. int operator>=(const real&x,const real&y)
  788. {return((y-x).sign()<=0?1:0);}
  789. int operator!=(const real&x,const real&y)
  790. {return((y-x).sign()!=0?1:0);}
  792. /*:46*/
  794. /*47:*/
  796. real fabs(real&x){ if (x.sign()<0) return -x; else return x; }
  797. real sq(const real&x){return x*x;}
  798. real hypot(const real&x,const real&y){return sqrt(x*x+y*y);}
  799. real powi(const real&x,int n)
  800. {real y= x,z= 1;
  801. int n_prefix= n;
  802. while(n_prefix>0)
  803. {
  804. if(n_prefix%2)z= z*y;
  805. n_prefix= n_prefix/2;
  806. y= y*y;
  807. }
  808. return z;
  809. }
  811. /*:47*/
  814. /*:12*/