;;;; -*- Mode: lisp -*-;;;;;;;; Copyright (c) 2007, 2008, 2011 Raymond Toy;;;;;;;; Permission is hereby granted, free of charge, to any person;;;; obtaining a copy of this software and associated documentation;;;; files (the "Software"), to deal in the Software without;;;; restriction, including without limitation the rights to use,;;;; copy, modify, merge, publish, distribute, sublicense, and/or sell;;;; copies of the Software, and to permit persons to whom the;;;; Software is furnished to do so, subject to the following;;;; conditions:;;;;;;;; The above copyright notice and this permission notice shall be;;;; included in all copies or substantial portions of the Software.;;;;;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES;;;; OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT;;;; HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,;;;; WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING;;;; FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR;;;; OTHER DEALINGS IN THE SOFTWARE.(in-package#:oct)(defmethodmake-qd((xcl:rational));; We should do something better than this.(make-instance'qd-real:value(rational-to-qdx)));; Determine which of x and y has the higher precision and return the;; value of the higher precision number. If both x and y are;; rationals, just return 1f0, for a single-float value.(defunfloat-contagion-2(xy)(etypecasex(cl:rational(etypecasey(cl:rational1f0)(cl:floaty)(qd-realy)))(single-float(etypecasey((orcl:rationalsingle-float)x)((ordouble-floatqd-real)y)))(double-float(etypecasey((orcl:rationalsingle-floatdouble-float)x)(qd-realy)))(qd-realx)));; Return a floating point (or complex) type of the highest precision;; among all of the given arguments.(defunfloat-contagion(&restargs);; It would be easy if we could just add the args together and let;; normal contagion do the work, but we could easily introduce;; overflows or other errors that way. So look at each argument and;; determine the precision and choose the highest precision. (etypecase(reduce#'float-contagion-2(mapcar#'realpart(if(cdrargs)args(list(carargs)0))))(single-float'single-float)(double-float'double-float)(qd-real'qd-real)))(defunapply-contagion(numberprecision)(etypecasenumber((orcl:realqd-real)(coercenumberprecision))((orcl:complexqd-complex)(complex(coerce(realpartnumber)precision)(coerce(imagpartnumber)precision)))));; WITH-FLOATING-POINT-CONTAGION - macro;;;; Determines the highest precision of the variables in VARLIST and;; converts each of the values to that precision.(defmacrowith-floating-point-contagion(varlist&bodybody)(let((precision(gensym"PRECISION-")))`(let((,precision(float-contagion,@varlist)))(let(,@(mapcar#'(lambda(v)`(,v(apply-contagion,v,precision)))varlist)),@body))))(defmethodadd1((anumber))(cl::1+a))(defmethodadd1((aqd-real))(make-instance'qd-real:value(add-qd-d(qd-valuea)1d0)))(defmethodsub1((anumber))(cl::1-a))(defmethodsub1((aqd-real))(make-instance'qd-real:value(sub-qd-d(qd-valuea)1d0)))(declaim(inline1+1-))(defun1+(x)(add1x))(defun1-(x)(sub1x))(defmethodtwo-arg-+((aqd-real)(bqd-real))(make-instance'qd-real:value(add-qd(qd-valuea)(qd-valueb))))(defmethodtwo-arg-+((aqd-real)(bcl:float))(make-instance'qd-real:value(add-qd-d(qd-valuea)(cl:floatb1d0))))#+cmu(defmethodtwo-arg-+((aqd-real)(bext:double-double-float))(make-instance'qd-real:value(add-qd-dd(qd-valuea)b)))(defmethodtwo-arg-+((areal)(bqd-real))(+ba))(defmethodtwo-arg-+((anumber)(bnumber))(cl:+ab))(defun+(&restargs)(if(nullargs)0(do((args(cdrargs)(cdrargs))(res(carargs)(two-arg-+res(carargs))))((nullargs)res))))(defmethodtwo-arg--((aqd-real)(bqd-real))(make-instance'qd-real:value(sub-qd(qd-valuea)(qd-valueb))))(defmethodtwo-arg--((aqd-real)(bcl:float))(make-instance'qd-real:value(sub-qd-d(qd-valuea)(cl:floatb1d0))))#+cmu(defmethodtwo-arg--((aqd-real)(bext:double-double-float))(make-instance'qd-real:value(sub-qd-dd(qd-valuea)b)))(defmethodtwo-arg--((acl:float)(bqd-real))(make-instance'qd-real:value(sub-d-qd(cl:floata1d0)(qd-valueb))))(defmethodtwo-arg--((anumber)(bnumber))(cl:-ab))(defmethodunary-minus((anumber))(cl:-a))(defmethodunary-minus((aqd-real))(make-instance'qd-real:value(neg-qd(qd-valuea))))(defun-(number&restmore-numbers)(ifmore-numbers(do((nlistmore-numbers(cdrnlist))(resultnumber))((atomnlist)result)(declare(listnlist))(setqresult(two-arg--result(carnlist))))(unary-minusnumber)))(defmethodtwo-arg-*((aqd-real)(bqd-real))(make-instance'qd-real:value(mul-qd(qd-valuea)(qd-valueb))))(defmethodtwo-arg-*((aqd-real)(bcl:float))(make-instance'qd-real:value(mul-qd-d(qd-valuea)(cl:floatb1d0))))#+cmu(defmethodtwo-arg-*((aqd-real)(bext:double-double-float));; We'd normally want to use mul-qd-dd, but mul-qd-dd is broken.(make-instance'qd-real:value(mul-qd(qd-valuea)(make-qd-ddb0w0))))(defmethodtwo-arg-*((areal)(bqd-real))(*ba))(defmethodtwo-arg-*((anumber)(bnumber))(cl:*ab))(defun*(&restargs)(if(nullargs)1(do((args(cdrargs)(cdrargs))(res(carargs)(two-arg-*res(carargs))))((nullargs)res))))(defmethodtwo-arg-/((aqd-real)(bqd-real))(make-instance'qd-real:value(div-qd(qd-valuea)(qd-valueb))))(defmethodtwo-arg-/((aqd-real)(bcl:float))(make-instance'qd-real:value(div-qd-d(qd-valuea)(cl:floatb1d0))))#+cmu(defmethodtwo-arg-/((aqd-real)(bext:double-double-float))(make-instance'qd-real:value(div-qd-dd(qd-valuea)b)))(defmethodtwo-arg-/((acl:float)(bqd-real))(make-instance'qd-real:value(div-qd(make-qd-d(cl:floata1d0))(qd-valueb))))#+cmu(defmethodtwo-arg-/((aext:double-double-float)(bqd-real))(make-instance'qd-real:value(div-qd(make-qd-dda0w0)(qd-valueb))))(defmethodtwo-arg-/((anumber)(bnumber))(cl:/ab))(defmethodunary-divide((anumber))(cl:/a))(defmethodunary-divide((aqd-real))(make-instance'qd-real:value(div-qd+qd-one+(qd-valuea))))(defun/(number&restmore-numbers)(ifmore-numbers(do((nlistmore-numbers(cdrnlist))(resultnumber))((atomnlist)result)(declare(listnlist))(setqresult(two-arg-/result(carnlist))))(unary-dividenumber)))(macrolet((frob(name&optional(type'real))(let((method-name(intern(concatenate'string(string'#:q)(symbol-namename))))(cl-name(intern(symbol-namename):cl))(qd-name(intern(concatenate'string(symbol-namename)(string'#:-qd)))))`(progn(defmethod,method-name((x,type))(,cl-namex))(defmethod,method-name((xqd-real))(,qd-name(qd-valuex)))(declaim(inline,name))(defun,name(x)(,method-namex))))))(frobzeropnumber)(frobplusp)(frobminusp))(defunbignum-to-qd(bignum)(make-instance'qd-real:value(rational-to-qdbignum)))(defunfloatp(x)(typepx'(orshort-floatsingle-floatdouble-floatlong-floatqd-real)))(defmethodqfloat((xreal)(num-typecl:float))(cl:floatxnum-type))(defmethodqfloat((xcl:float)(num-typeqd-real))(make-instance'qd-real:value(make-qd-d(cl:floatx1d0))))(defmethodqfloat((xinteger)(num-typeqd-real))(cond((typepx'fixnum)(make-instance'qd-real:value(make-qd-d(cl:floatx1d0))))(t;; A bignum(bignum-to-qdx))))#+nil(defmethodqfloat((xratio)(num-typeqd-real));; This probably has some issues with roundoff(two-arg-/(qfloat(numeratorx)num-type)(qfloat(denominatorx)num-type)))(defmethodqfloat((xratio)(num-typeqd-real))(make-instance'qd-real:value(rational-to-qdx)))#+cmu(defmethodqfloat((xext:double-double-float)(num-typeqd-real))(make-instance'qd-real:value(make-qd-ddx0w0)))(defmethodqfloat((xqd-real)(num-typecl:float))(multiple-value-bind(q0q1q2q3)(qd-parts(qd-valuex))(cl:float(cl:+q0q1q2q3)num-type)))#+cmu(defmethodqfloat((xqd-real)(num-typeext:double-double-float))(multiple-value-bind(q0q1q2q3)(qd-parts(qd-valuex))(cl:+(cl:floatq01w0)(cl:floatq11w0)(cl:floatq21w0)(cl:floatq31w0))))(defmethodqfloat((xqd-real)(num-typeqd-real))x)(declaim(inlinefloat))(defunfloat(x&optional(othernilotherp))(ifotherp(qfloatxother)(if(floatpx)x(qfloatx1.0))))(defmethodqrealpart((xnumber))(cl:realpartx))(defmethodqrealpart((xqd-real))x)(defmethodqrealpart((xqd-complex))(make-instance'qd-real:value(qd-realx)))(defunrealpart(x)(qrealpartx))(defmethodqimagpart((xnumber))(cl:imagpartx))(defmethodqimagpart((xqd-real))(make-qd0d0))(defmethodqimagpart((xqd-complex))(make-instance'qd-real:value(qd-imagx)))(defunimagpart(x)(qimagpartx))(defmethodqconjugate((anumber))(cl:conjugatea))(defmethodqconjugate((aqd-real))(make-instance'qd-real:value(qd-valuea)))(defmethodqconjugate((aqd-complex))(make-instance'qd-complex:real(qd-reala):imag(neg-qd(qd-imaga))))(defunconjugate(z)(qconjugatez))(defmethodqscale-float((fcl:float)(ninteger))(cl:scale-floatfn))(defmethodqscale-float((fqd-real)(ninteger))(make-instance'qd-real:value(scale-float-qd(qd-valuef)n)))(declaim(inlinescale-float))(defunscale-float(fn)(qscale-floatfn))(macrolet((frob(op)(let((method(intern(concatenate'string(string'#:two-arg-)(symbol-nameop))))(cl-fun(find-symbol(symbol-nameop):cl))(qd-fun(intern(concatenate'string(string'#:qd-)(symbol-nameop))'#:octi)))`(progn(defmethod,method((areal)(breal))(,cl-funab))(defmethod,method((aqd-real)(breal))(,qd-fun(qd-valuea)(make-qd-d(cl:floatb1d0))))(defmethod,method((areal)(bqd-real));; This is not really right if A is a rational. We're;; supposed to compare them as rationals.(,qd-fun(make-qd-d(cl:floata1d0))(qd-valueb)))(defmethod,method((aqd-real)(bqd-real))(,qd-fun(qd-valuea)(qd-valueb)))(defun,op(number&restmore-numbers)"Returns T if its arguments are in strictly increasing order, NIL otherwise."(declare(optimize(safety2))(dynamic-extentmore-numbers))(do*((nnumber(carnlist))(nlistmore-numbers(cdrnlist)))((atomnlist)t)(declare(listnlist))(if(not(,methodn(carnlist)))(returnnil))))))))(frob<)(frob>)(frob<=)(frob>=));; Handle the special functions for a real argument. Complex args are;; handled elsewhere.(macrolet((frob(name)(let((method-name(intern(concatenate'string(string'#:q)(symbol-namename))))(cl-name(intern(symbol-namename):cl))(qd-name(intern(concatenate'string(symbol-namename)(string'#:-qd)))))`(progn(defmethod,name((xnumber))(,cl-namex))(defmethod,name((xqd-real))(make-instance'qd-real:value(,qd-name(qd-valuex))))))))(frobabs)(frobexp)(frobsin)(frobcos)(frobtan);;(frob asin);;(frob acos)(frobsinh)(frobcosh)(frobtanh)(frobasinh);;(frob acosh);;(frob atanh))(defmethodsqrt((xnumber))(cl:sqrtx))(defmethodsqrt((xqd-real))(if(minuspx)(make-instance'qd-complex:real+qd-zero+:imag(sqrt-qd(neg-qd(qd-valuex))))(make-instance'qd-real:value(sqrt-qd(qd-valuex)))))(defunscalb(xn)"Compute 2^N * X without compute 2^N first (use properties of theunderlying floating-point format"(declare(typeqd-realx))(scale-floatxn))(declaim(inlineqd-cssqs))(defunqd-cssqs(z)(multiple-value-bind(rhok)(octi::hypot-aux-qd(qd-value(realpartz))(qd-value(imagpartz)))(values(make-instance'qd-real:valuerho)k)))#+nil(defmethodqabs((zqd-complex));; sqrt(x^2+y^2);; If |x| > |y| then sqrt(x^2+y^2) = |x|*sqrt(1+(y/x)^2)(multiple-value-bind(abs^2rho)(hypot-qd(qd-value(realpartz))(qd-value(imagpartz)))(scale-float(make-instance'qd-real:value(sqrtabs^2))rho)))(defmethodabs((zqd-complex));; sqrt(x^2+y^2);; If |x| > |y| then sqrt(x^2+y^2) = |x|*sqrt(1+(y/x)^2)(make-instance'qd-real:value(hypot-qd(qd-value(realpartz))(qd-value(imagpartz)))))(defmethodlog((anumber)&optionalb)(ifb(cl:logab)(cl:loga)))(defmethodlog((aqd-real)&optionalb)(ifb(/(loga)(logb))(if(minusp(float-signa))(make-instance'qd-complex:real(log-qd(abs-qd(qd-valuea))):imag+qd-pi+)(make-instance'qd-real:value(log-qd(qd-valuea))))))(defmethodlog1p((aqd-real))(make-instance'qd-real:value(log1p-qd(qd-valuea))))(defmethodatan((yreal)&optionalx)(cond(x(cond((typepx'qd-real)(make-instance'qd-real:value(atan2-qd(qd-valuey)(qd-valuex))))(t(cl:atanyx))))(t(cl:atany))))(defmethodatan((yqd-real)&optionalx)(make-instance'qd-real:value(ifx(atan2-qd(qd-valuey)(qd-valuex))(atan-qd(qd-valuey)))))(defmethodqexpt((xnumber)(ynumber))(cl:exptxy))(defmethodqexpt((xnumber)(yqd-real))(exp(*y(log(apply-contagionx'qd-real)))))(defmethodqexpt((xnumber)(yqd-complex))(exp(*y(log(apply-contagionx'qd-real)))))(defmethodqexpt((xqd-real)(yreal))(exp(*y(logx))))(defmethodqexpt((xqd-real)(ycl:complex))(exp(*(make-instance'qd-complex:real(qd-value(realparty)):imag(qd-value(imagparty)))(logx))))(defmethodqexpt((xqd-real)(yqd-real));; x^y = exp(y*log(x))(exp(*y(logx))))(defmethodqexpt((xqd-real)(yinteger))(make-instance'qd-real:value(npow(qd-valuex)y)))(declaim(inlineexpt))(defunexpt(xy)(qexptxy))(defmethodtwo-arg-=((anumber)(bnumber))(cl:=ab))(defmethodtwo-arg-=((aqd-real)(bnumber))(if(cl:realpb)(qd-=(qd-valuea)(make-qd-d(cl:floatb1d0)))nil))(defmethodtwo-arg-=((anumber)(bqd-real))(if(cl:realpa)(qd-=(make-qd-d(cl:floata1d0))(qd-valueb))nil))(defmethodtwo-arg-=((aqd-complex)b)(and(two-arg-=(realparta)(realpartb))(two-arg-=(imagparta)(imagpartb))))(defmethodtwo-arg-=(a(bqd-complex))(and(two-arg-=(realparta)(realpartb))(two-arg-=(imagparta)(imagpartb))))(defmethodtwo-arg-=((aqd-real)(bqd-real))(qd-=(qd-valuea)(qd-valueb)))(defun=(number&restmore-numbers)"Returns T if all of its arguments are numerically equal, NIL otherwise."(declare(optimize(safety2))(dynamic-extentmore-numbers))(do((nlistmore-numbers(cdrnlist)))((atomnlist)t)(declare(listnlist))(if(not(two-arg-=(carnlist)number))(returnnil))))(defun/=(number&restmore-numbers)"Returns T if no two of its arguments are numerically equal, NIL otherwise."(declare(optimize(safety2))(dynamic-extentmore-numbers))(do*((headnumber(carnlist))(nlistmore-numbers(cdrnlist)))((atomnlist)t)(declare(listnlist))(unless(do*((nlnlist(cdrnl)))((atomnl)t)(declare(listnl))(if(two-arg-=head(carnl))(returnnil)))(returnnil))))(defmethodqcomplex((xcl:real)(ycl:real))(cl:complexxy))(defmethodqcomplex((xcl:real)(yqd-real))(qcomplex(make-qdx)y))(defmethodqcomplex((xqd-real)(yqd-real))(make-instance'qd-complex:real(qd-valuex):imag(qd-valuey)))(defmethodqcomplex((xqd-real)(ycl:real))(make-instance'qd-complex:real(qd-valuex):imag(make-qd-dy)))(defuncomplex(x&optional(y0))(qcomplexxy))(defmethodqinteger-decode-float((fcl:float))(cl:integer-decode-floatf))(defmethodqinteger-decode-float((fqd-real))(integer-decode-qd(qd-valuef)))(declaim(inlineinteger-decode-float))(defuninteger-decode-float(f)(qinteger-decode-floatf))(defmethodqdecode-float((fcl:float))(cl:decode-floatf))(defmethodqdecode-float((fqd-real))(multiple-value-bind(fracexps)(decode-float-qd(qd-valuef))(values(make-instance'qd-real:valuefrac)exp(make-instance'qd-real:values))))(declaim(inlinedecode-float))(defundecode-float(f)(qdecode-floatf))(defmethodqfloor((xreal)&optionaly)(ify(cl:floorxy)(cl:floorx)))(defmethodqfloor((xqd-real)&optionaly)(if(andy(/=y1))(let((f(qfloor(/xy))))(valuesf(-x(*fy))))(let((f(ffloor-qd(qd-valuex))))(multiple-value-bind(intexpsign)(integer-decode-qdf)(values(ash(*signint)exp)(make-instance'qd-real:value(qd-value(-x(make-instance'qd-real:valuef)))))))))(defunfloor(x&optionaly)(qfloorxy))(defmethodqffloor((xreal)&optionaly)(ify(cl:ffloorxy)(cl:ffloorx)))(defmethodqffloor((xqd-real)&optionaly)(if(andy(/=y1))(let((f(qffloor(/xy))))(valuesf(-x(*fy))))(let((f(make-instance'qd-real:value(ffloor-qd(qd-valuex)))))(valuesf(-xf)))))(defunffloor(x&optionaly)(qffloorxy))(defunceiling(x&optionaly)(multiple-value-bind(frem)(floorxy)(if(zeroprem)(valuesfrem)(values(+f1)(-rem1)))))(defunfceiling(x&optionaly)(multiple-value-bind(frem)(ffloorxy)(if(zeroprem)(valuesfrem)(values(+f1)(-rem1)))))(defuntruncate(x&optional(y1))(if(minuspx)(ceilingxy)(floorxy)))(defunrem(xy)(nth-value1(truncatexy)))(defunmod(xy)(nth-value1(floorxy)))(defunftruncate(x&optional(y1))(if(minuspx)(fceilingxy)(ffloorxy)))(defmethod%unary-round((xreal))(cl::roundx))(defmethod%unary-round((numberqd-real))(multiple-value-bind(bitsexp)(integer-decode-floatnumber)(let*((shifted(ashbitsexp))(rounded(if(and(minuspexp)(oddpshifted)(not(zerop(logandbits(ash1(--1exp))))))(1+shifted)shifted)))(if(minuspnumber)(-rounded)rounded))))(defunround(number&optional(divisor1))(if(eqldivisor1)(let((r(%unary-roundnumber)))(valuesr(-numberr)))(multiple-value-bind(trurem)(truncatenumberdivisor)(if(zeroprem)(valuestrurem)(let((thresh(/(absdivisor)2)))(cond((or(>remthresh)(and(=remthresh)(oddptru)))(if(minuspdivisor)(values(-tru1)(+remdivisor))(values(+tru1)(-remdivisor))))((let((-thresh(-thresh)))(or(<rem-thresh)(and(=rem-thresh)(oddptru))))(if(minuspdivisor)(values(+tru1)(-remdivisor))(values(-tru1)(+remdivisor))))(t(valuestrurem))))))))(defunfround(number&optional(divisor1))"Same as ROUND, but returns first value as a float."(multiple-value-bind(resrem)(roundnumberdivisor)(values(floatres(if(floatprem)rem1.0))rem)))(defmethodqfloat-sign((areal)&optional(f(float1a)))(cl:float-signaf))(defmethodqfloat-sign((aqd-real)&optionalf)(iff(make-instance'qd-real:value(mul-qd-d(abs-qd(qd-valuef))(cl:float-sign(qd-0(qd-valuea)))))(make-instance'qd-real:value(make-qd-d(cl:float-sign(qd-0(qd-valuea)))))))(declaim(inlinefloat-sign))(defunfloat-sign(n&optional(float2nilfloat2p))(iffloat2p(qfloat-signnfloat2)(qfloat-signn)))(defunmax(number&restmore-numbers)"Returns the greatest of its arguments."(declare(optimize(safety2))(type(orrealqd-real)number)(dynamic-extentmore-numbers))(dolist(realmore-numbers)(when(>realnumber)(setqnumberreal)))number)(defunmin(number&restmore-numbers)"Returns the least of its arguments."(declare(optimize(safety2))(type(orrealqd-real)number)(dynamic-extentmore-numbers))(do((nlistmore-numbers(cdrnlist))(result(the(orrealqd-real)number)))((nullnlist)(returnresult))(declare(listnlist))(if(<(carnlist)result)(setqresult(carnlist)))))(defmethodasin((xnumber))(cl:asinx))(defmethodasin((xqd-real))(if(<=-1x1)(make-instance'qd-real:value(asin-qd(qd-valuex)))(qd-complex-asinx)))(defmethodacos((xnumber))(cl:acosx))(defmethodacos((xqd-real))(cond((>(absx)1)(qd-complex-acosx))(t(make-instance'qd-real:value(acos-qd(qd-valuex))))))(defmethodacosh((xnumber))(cl:acoshx))(defmethodacosh((xqd-real))(if(<x1)(qd-complex-acoshx)(make-instance'qd-real:value(acosh-qd(qd-valuex)))))(defmethodatanh((xnumber))(cl:atanhx))(defmethodatanh((xqd-real))(if(>(absx)1)(qd-complex-atanhx)(make-instance'qd-real:value(atanh-qd(qd-valuex)))))(defmethodcis((xreal))(cl:cisx))(defmethodcis((xqd-real))(multiple-value-bind(sc)(sincos-qd(qd-valuex))(make-instance'qd-complex:realc:imags)))(defmethodphase((xnumber))(cl:phasex))(defmethodphase((xqd-real))(if(minuspx)(-+pi+)(make-instance'qd-real:value(make-qd-d0d0))))(defunsignum(number)"If NUMBER is zero, return NUMBER, else return (/ NUMBER (ABS NUMBER))."(if(zeropnumber)number(if(rationalpnumber)(if(pluspnumber)1-1)(/number(absnumber)))))(defmethodrandom((xcl:real)&optional(state*random-state*))(cl:randomxstate))(defmethodrandom((xqd-real)&optional(state*random-state*))(*x(make-instance'qd-real:value(octi:random-qdstate))))(defmethodfloat-digits((xcl:real))(cl:float-digitsx))(defmethodfloat-digits((xqd-real))(*4(float-digits1d0)))(defmethodrational((xreal))(cl:rationalx))(defmethodrational((xqd-real))(with-qd-parts(x0x1x2x3)(qd-valuex)(+(cl:rationalx0)(cl:rationalx1)(cl:rationalx2)(cl:rationalx3))))(defmethodrationalize((xcl:real))(cl:rationalizex));;; The algorithm here is the method described in CLISP. Bruno Haible has;;; graciously given permission to use this algorithm. He says, "You can use;;; it, if you present the following explanation of the algorithm.";;;;;; Algorithm (recursively presented):;;; If x is a rational number, return x.;;; If x = 0.0, return 0.;;; If x < 0.0, return (- (rationalize (- x))).;;; If x > 0.0:;;; Call (integer-decode-float x). It returns a m,e,s=1 (mantissa,;;; exponent, sign).;;; If m = 0 or e >= 0: return x = m*2^e.;;; Search a rational number between a = (m-1/2)*2^e and b = (m+1/2)*2^e;;; with smallest possible numerator and denominator.;;; Note 1: If m is a power of 2, we ought to take a = (m-1/4)*2^e.;;; But in this case the result will be x itself anyway, regardless of;;; the choice of a. Therefore we can simply ignore this case.;;; Note 2: At first, we need to consider the closed interval [a,b].;;; but since a and b have the denominator 2^(|e|+1) whereas x itself;;; has a denominator <= 2^|e|, we can restrict the seach to the open;;; interval (a,b).;;; So, for given a and b (0 < a < b) we are searching a rational number;;; y with a <= y <= b.;;; Recursive algorithm fraction_between(a,b):;;; c := (ceiling a);;; if c < b;;; then return c ; because a <= c < b, c integer;;; else;;; ; a is not integer (otherwise we would have had c = a < b);;; k := c-1 ; k = floor(a), k < a < b <= k+1;;; return y = k + 1/fraction_between(1/(b-k), 1/(a-k));;; ; note 1 <= 1/(b-k) < 1/(a-k);;;;;; You can see that we are actually computing a continued fraction expansion.;;;;;; Algorithm (iterative):;;; If x is rational, return x.;;; Call (integer-decode-float x). It returns a m,e,s (mantissa,;;; exponent, sign).;;; If m = 0 or e >= 0, return m*2^e*s. (This includes the case x = 0.0.);;; Create rational numbers a := (2*m-1)*2^(e-1) and b := (2*m+1)*2^(e-1);;; (positive and already in lowest terms because the denominator is a;;; power of two and the numerator is odd).;;; Start a continued fraction expansion;;; p[-1] := 0, p[0] := 1, q[-1] := 1, q[0] := 0, i := 0.;;; Loop;;; c := (ceiling a);;; if c >= b;;; then k := c-1, partial_quotient(k), (a,b) := (1/(b-k),1/(a-k)),;;; goto Loop;;; finally partial_quotient(c).;;; Here partial_quotient(c) denotes the iteration;;; i := i+1, p[i] := c*p[i-1]+p[i-2], q[i] := c*q[i-1]+q[i-2].;;; At the end, return s * (p[i]/q[i]).;;; This rational number is already in lowest terms because;;; p[i]*q[i-1]-p[i-1]*q[i] = (-1)^i.;;;(defmethodrationalize((xqd-real));; This is a fairly straigtforward implementation of the iterative;; algorithm above.(multiple-value-bind(fracexposign)(integer-decode-floatx)(cond((or(zeropfrac)(>=expo0))(if(minuspsign)(-(ashfracexpo))(ashfracexpo)))(t;; expo < 0 and (2*m-1) and (2*m+1) are coprime to 2^(1-e),;; so build the fraction up immediately, without having to do;; a gcd.(let((a(/(-(*2frac)1)(ash1(-1expo))))(b(/(+(*2frac)1)(ash1(-1expo))))(p00)(q01)(p11)(q10))(do((c(ceilinga)(ceilinga)))((<cb)(let((top(+(*cp1)p0))(bot(+(*cq1)q0)))(/(if(minuspsign)(-top)top)bot)))(let*((k(-c1))(p2(+(*kp1)p0))(q2(+(*kq1)q0)))(psetfa(/(-bk))b(/(-ak)))(setfp0p1q0q1p1p2q1q2))))))))(define-compiler-macro+(&wholeform&restargs)(declare(ignoreform))(if(nullargs)0(do((args(cdrargs)(cdrargs))(res(carargs)`(two-arg-+,res,(carargs))))((nullargs)res))))(define-compiler-macro-(&wholeformnumber&restmore-numbers)(declare(ignoreform))(ifmore-numbers(do((nlistmore-numbers(cdrnlist))(resultnumber))((atomnlist)result)(declare(listnlist))(setqresult`(two-arg--,result,(carnlist))))`(unary-minus,number)))(define-compiler-macro*(&wholeform&restargs)(declare(ignoreform))(if(nullargs)1(do((args(cdrargs)(cdrargs))(res(carargs)`(two-arg-*,res,(carargs))))((nullargs)res))))(define-compiler-macro/(number&restmore-numbers)(ifmore-numbers(do((nlistmore-numbers(cdrnlist))(resultnumber))((atomnlist)result)(declare(listnlist))(setqresult`(two-arg-/,result,(carnlist))))`(unary-divide,number)));; Compiler macros to convert <, >, <=, and >= into multiple calls of;; the corresponding two-arg-<foo> function.(macrolet((frob(op)(let((method(intern(concatenate'string(string'#:two-arg-)(symbol-nameop)))))`(define-compiler-macro,op(number&restmore-numbers)(do*((nnumber(carnlist))(nlistmore-numbers(cdrnlist))(resnil))((atomnlist)`(and,@(nreverseres)))(push`(,',method,n,(carnlist))res))))))(frob<)(frob>)(frob<=)(frob>=))(define-compiler-macro/=(&wholeformnumber&restmore-numbers);; Convert (/= x y) to (not (two-arg-= x y)). Should we try to;; handle a few more cases?(if(cdrmore-numbers)form`(not(two-arg-=,number,(carmore-numbers)))));; Define compiler macro the convert two-arg-foo into the appropriate;; CL function or QD-REAL function so we don't have to do CLOS;; dispatch.#+(or)(macrolet((frob(namecl-opqd-op)`(define-compiler-macro,name(&wholeformxy&environmentenv)(flet((arg-type(arg)(multiple-value-bind(def-typelocalpdecl)(ext:variable-informationargenv)(declare(ignorelocalp))(whendef-type(cdr(assoc'typedecl))))))(let((x-type(arg-typex))(y-type(arg-typey)))(cond((and(subtypepx-type'cl:number)(subtypepy-type'cl:number))`(,',cl-op,x,y))((and(subtypepx-type'qd-real)(subtypepy-type'qd-real))`(make-instance'qd-real:value(,',qd-op(qd-value,x)(qd-value,y))))(t;; Don't know how to handle this, so give up.form)))))))(frobtwo-arg-+cl:+add-qd)(frobtwo-arg--cl:-sub-qd)(frobtwo-arg-*cl:*mul-qd)(frobtwo-arg-/cl:/div-qd))#+(or)(macrolet((frob(namecl-opqd-opcl-qd-opqd-cl-op)`(define-compiler-macro,name(&wholeformxy&environmentenv)(flet((arg-type(arg)(multiple-value-bind(def-typelocalpdecl)(ext:variable-informationargenv)(declare(ignorelocalp))(whendef-type(cdr(assoc'typedecl))))))(let((x-type(arg-typex))(y-type(arg-typey)))(cond((subtypepx-type'cl:float)(cond((subtypepy-type'cl:number)`(,',cl-op,x,y))((subtypepy-type'qd-real)(if,cl-qd-op`(make-instance'qd-real:value(,',cl-qd-op(cl:float,x1d0)(qd-value,y)))form))(tform)))((subtypepx-type'qd-real)(cond((subtypepy-type'cl:float)(if,qd-cl-op`(make-instance'qd-real:value(,',qd-cl-op(qd-value,x)(float,y1d0)))form))((subtypepy-type'qd-real)`(make-instance'qd-real:value(,',qd-op(qd-value,x)(qd-value,y))))(tform)))(t;; Don't know how to handle this, so give up.form)))))))(frobtwo-arg-+cl:+add-qdadd-d-qdadd-qd-d)(frobtwo-arg--cl:-sub-qdsub-d-qdsub-qd-d)(frobtwo-arg-*cl:*mul-qdmul-d-qdmul-qd-d)(frobtwo-arg-/cl:/div-qdnilnil))(defgenericepsilon(m)(:documentation"Return an epsilon value of the same precision as the argument. It isthe smallest number x such that 1+x /= x. The argument can becomplex"))(defmethodepsilon((mcl:float))(etypecasem(single-floatsingle-float-epsilon)(double-floatdouble-float-epsilon)))(defmethodepsilon((mcl:complex))(epsilon(realpartm)))(defmethodepsilon((mqd-real));; What is the epsilon value for a quad-double? This is complicated;; by the fact that things like (+ #q1 #q1q-100) is representable as;; a quad-double. For most purposes we want epsilon to be close to;; the 212 bits of precision (4*53 bits) that we normally have with;; a quad-double.(scale-float+qd-real-one+-212))(defmethodepsilon((mqd-complex))(epsilon(realpartm)))(defgenericfloat-pi(x)(:documentation"Return a floating-point value of the mathematical constant pi that isthe same precision as the argument. The argument can be complex."))(defmethodfloat-pi((xcl:rational))(floatpi1f0))(defmethodfloat-pi((xcl:float))(floatpix))(defmethodfloat-pi((xqd-real))+pi+)(defmethodfloat-pi((zcl:complex))(floatpi(realpartz)))(defmethodfloat-pi((zqd-complex))+pi+)(defmethodfloat-nan-p((xcl:float));; CMUCL has ext:float-nan-p. Should we use that instead?(not(=xx)))(defmethodfloat-nan-p((xqd-real))(float-nan-p(qd-parts(qd-valuex))))(define-conditiondomain-error(simple-error)((function-name:accessorcondition-function-name:initarg:function-name))(:report(lambda(conditionstream)(formatstream"Domain Error for function ~S:~&"(condition-function-namecondition))(pprint-logical-block(streamnil:per-line-prefix" ")(apply#'formatstream(simple-condition-format-controlcondition)(simple-condition-format-argumentscondition))))))