![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to math.sml CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [Research Unix] / researchv10no / cmd / sml / src / boot|
Return to math.sml CVS log | Up to [Research Unix] / researchv10no / cmd / sml / src / boot |
1.1 ! root 1: (* Copyright 1989 by AT&T Bell Laboratories *) ! 2: (* The following functions were adapted from the 4.3BSD math library. ! 3: Eventually, each machine supported should have a hand-coded math ! 4: functor with more efficient versions of these functions. ! 5: ! 6: *************************************************************************** ! 7: * * ! 8: * Copyright (c) 1985 Regents of the University of California. * ! 9: * * ! 10: * Use and reproduction of this software are granted in accordance with * ! 11: * the terms and conditions specified in the Berkeley Software License * ! 12: * Agreement (in particular, this entails acknowledgement of the programs' * ! 13: * source, and inclusion of this notice) with the additional understanding * ! 14: * that all recipients should regard themselves as participants in an * ! 15: * ongoing research project and hence should feel obligated to report * ! 16: * their experiences (good or bad) with these elementary function codes, * ! 17: * using "sendbug 4bsd-bugs@BERKELEY", to the authors. * ! 18: * * ! 19: * K.C. Ng, with Z-S. Alex Liu, S. McDonald, P. Tang, W. Kahan. * ! 20: * Revised on 5/10/85, 5/13/85, 6/14/85, 8/20/85, 8/27/85, 9/11/85. * ! 21: * * ! 22: *************************************************************************** ! 23: ! 24: *) ! 25: ! 26: signature MATH = ! 27: sig ! 28: exception Sqrt and Exp and Ln ! 29: val sqrt : real -> real ! 30: val sin : real -> real ! 31: val cos : real -> real ! 32: val arctan : real -> real ! 33: val exp : real -> real ! 34: val ln : real -> real ! 35: end ! 36: ! 37: structure Math : MATH = struct ! 38: ! 39: structure Assembly : ASSEMBLY = Core.Assembly ! 40: ! 41: infix 7 * / ! 42: infix 6 + - ! 43: infix 4 > < >= <= ! 44: ! 45: exception Real = Assembly.Real ! 46: exception Sqrt ! 47: exception Exp ! 48: exception Ln ! 49: ! 50: val scalb = Assembly.A.scalb ! 51: val logb = Assembly.A.logb ! 52: ! 53: val ~ = InLine.fneg ! 54: val op + = InLine.fadd ! 55: val op - = InLine.fsub ! 56: val op * = InLine.fmul ! 57: val op / = InLine.fdiv ! 58: val op > = InLine.fgt ! 59: val op < = InLine.flt ! 60: val op >= = InLine.fge ! 61: val op <= = InLine.fle ! 62: ! 63: val ieql = InLine.ieql ! 64: val ineq = InLine.ineq ! 65: val feql = InLine.feql ! 66: ! 67: val negone = ~1.0 ! 68: val zero = 0.0 ! 69: val half = 0.5 ! 70: val one = 1.0 ! 71: val two = 2.0 ! 72: val four = 4.0 ! 73: ! 74: fun copysign(a,b) = ! 75: case (a<zero,b<zero) of ! 76: (true,true) => a ! 77: | (false,false) => a ! 78: | _ => ~a ! 79: ! 80: fun abs x = if x < zero then ~x else x ! 81: fun mod(a,b) = InLine.-(a,InLine.*(InLine.div(a,b),b)) ! 82: local fun rtoi (pred,init) = ! 83: let fun loop n = ! 84: if pred n ! 85: then init ! 86: else let val (i,r) = loop(n * half) ! 87: val i' = InLine.lshift(i,1) ! 88: val r' = r+r ! 89: in if n-r' < one then (i',r') else (InLine.+(i',1),r'+one) ! 90: end ! 91: in loop ! 92: end ! 93: val pton = rtoi(fn x => x < one, (0,zero)) ! 94: val nton = rtoi(fn x => x >= negone, (~1,negone)) ! 95: fun fl n = if n < zero then nton n else pton n ! 96: fun tr n = #1(pton n) ! 97: in ! 98: fun truncate n = if n < zero then InLine.~(tr(~n)) else tr n ! 99: val floor = Assembly.A.floor ! 100: fun realfloor n = #2(fl n) ! 101: fun ceiling n = InLine.~(floor(~n)) ! 102: end ! 103: local fun loop 0 = zero ! 104: | loop n = let val x = two * loop(InLine.rshift(n,1)) ! 105: in case InLine.andb(n,1) of 0 => x | 1 => x + one ! 106: end ! 107: in fun real n = if InLine.<(n,0) then ~(loop(InLine.~ n)) else loop n ! 108: end ! 109: ! 110: (* Not a true drem - rounding on .5 should go to the even case. *) ! 111: fun drem(x,y) = ! 112: let val n = x/y+half ! 113: val N = if abs n < 9.007199254741E15 (* about 2^53 *) ! 114: then realfloor n ! 115: handle Real _ => ! 116: scalb(realfloor(scalb(n,InLine.~(logb n))),logb n) ! 117: else n ! 118: in x - N * y ! 119: end ! 120: ! 121: (* sin/cos *) ! 122: local ! 123: val S0 = ~1.6666666666666463126E~1 ! 124: val S1 = 8.3333333332992771264E~3 ! 125: val S2 = ~1.9841269816180999116E~4 ! 126: val S3 = 2.7557309793219876880E~6 ! 127: val S4 = ~2.5050225177523807003E~8 ! 128: val S5 = 1.5868926979889205164E~10 ! 129: in fun sin__S z = (z*(S0+z*(S1+z*(S2+z*(S3+z*(S4+z*S5)))))) ! 130: end ! 131: ! 132: local ! 133: val C0 = 4.1666666666666504759E~2 ! 134: val C1 = ~1.3888888888865301516E~3 ! 135: val C2 = 2.4801587269650015769E~5 ! 136: val C3 = ~2.7557304623183959811E~7 ! 137: val C4 = 2.0873958177697780076E~9 ! 138: val C5 = ~1.1250289076471311557E~11 ! 139: in fun cos__C z = (z*z*(C0+z*(C1+z*(C2+z*(C3+z*(C4+z*C5)))))) ! 140: end ! 141: ! 142: val PIo4 = 7.8539816339744827900E~1 ! 143: val PIo2 = 1.5707963267948965580E0 ! 144: val PI3o4 = 2.3561944901923448370E0 ! 145: val PI = 3.1415926535897931160E0 ! 146: val PI2 = 6.2831853071795862320E0 ! 147: ! 148: local ! 149: val thresh = 2.6117239648121182150E~1 ! 150: in fun S y = y + y * sin__S(y*y) ! 151: fun C y = ! 152: let val yy = y*y ! 153: val c = cos__C yy ! 154: val Y = yy/two ! 155: in if Y < thresh then one - (Y - c) ! 156: else half - (Y - half - c) ! 157: end ! 158: end ! 159: fun sin x = ! 160: let val x = drem(x,PI2) ! 161: val a = abs x ! 162: val y = if a >= PI3o4 then copysign(PI-a,x) else a ! 163: in if y >= PIo4 then copysign(C(PIo2-a),x) ! 164: else S y ! 165: end ! 166: ! 167: fun cos x = ! 168: let val x = drem(x,PI2) ! 169: val a = abs x ! 170: val (s,y) = if a >= PI3o4 ! 171: then (fn x => ~x,PI-a) ! 172: else (fn x => x,a) ! 173: in if y < PIo4 then s(C y) ! 174: else S(PIo2-a) ! 175: end ! 176: ! 177: local ! 178: val p1 = 1.3887401997267371720E~2 ! 179: val p2 = 3.3044019718331897649E~5 ! 180: val q1 = 1.1110813732786649355E~1 ! 181: val q2 = 9.9176615021572857300E~4 ! 182: in fun exp__E(x:real,c:real) = ! 183: let val small=1.0E~19 ! 184: val x = abs x ! 185: val z = x*x ! 186: val p = z*(p1+z*p2) ! 187: val q = z*(q1+z*q2) ! 188: val xp= x*p ! 189: val xh= x*half ! 190: val w = xh-(q-xp) ! 191: val p = p+p ! 192: val c = c+x*((xh*w-(q-(p+xp)))/(one-w)+c) ! 193: in if(copysign(x,one)>small) then z*half+c ! 194: else copysign(zero,x) (* could be inexact *) ! 195: end ! 196: end ! 197: ! 198: (* for exp and ln *) ! 199: val ln2hi = 6.9314718036912381649E~1 ! 200: val ln2lo = 1.9082149292705877000E~10 ! 201: val sqrt2 = 1.4142135623730951455E0 ! 202: val lnhuge = 7.1602103751842355450E2 ! 203: val lntiny = ~7.5137154372698068983E2 ! 204: val invln2 = 1.4426950408889633870E0 ! 205: ! 206: fun exp(x:real) = ! 207: let fun exp_norm() = ! 208: let (* argument reduction : x --> x - k*ln2 *) ! 209: val k = floor(invln2*x+copysign(half,x)) (* k=NINT(x/ln2) *) ! 210: val K = real k ! 211: (* express x-k*ln2 as z+c *) ! 212: val hi = x-K*ln2hi ! 213: val lo = K*ln2lo ! 214: val z = hi - lo ! 215: val c = (hi-z)-lo ! 216: (* return 2^k*[expm1(x) + 1] *) ! 217: val z = z + exp__E(z,c) ! 218: in scalb(z+one,k) ! 219: end ! 220: in if x <= lnhuge then if x >= lntiny ! 221: then exp_norm() handle Real _ => raise Exp ! 222: else zero (* exp(-big#) underflows to zero *) ! 223: (* exp(-INF) is zero *) ! 224: else raise Exp (* exp(INF) is INF, exp(+big#) overflows to INF *) ! 225: end ! 226: ! 227: local ! 228: val L1 = 6.6666666666667340202E~1 ! 229: val L2 = 3.9999999999416702146E~1 ! 230: val L3 = 2.8571428742008753154E~1 ! 231: val L4 = 2.2222198607186277597E~1 ! 232: val L5 = 1.8183562745289935658E~1 ! 233: val L6 = 1.5314087275331442206E~1 ! 234: val L7 = 1.4795612545334174692E~1 ! 235: in fun log__L(z) = z*(L1+z*(L2+z*(L3+z*(L4+z*(L5+z*(L6+z*L7)))))) ! 236: end ! 237: ! 238: fun ln(x:real) = ! 239: let fun findln() = ! 240: let (* argument reduction *) ! 241: val k = logb(x) ! 242: val x = scalb(x,InLine.~ k) ! 243: val (x,k) = if ieql(k,~1022) (* subnormal no. *) ! 244: then let val n = logb(x) ! 245: in (scalb(x,InLine.~ n),InLine.+(k,n)) end ! 246: else (x,k) ! 247: val (k,x) = if x >= sqrt2 then (InLine.+(k,1),x*half) else (k,x) ! 248: val K = real k ! 249: val x = x - one ! 250: (* compute log(1+x) *) ! 251: val s = x/(two+x) ! 252: val t = x*x*half ! 253: val z = K*ln2lo+s*(t+log__L(s*s)) ! 254: val x = x + (z - t) ! 255: in K*ln2hi+x end ! 256: in if x <= zero then raise Ln ! 257: else findln() handle Real _ => raise Ln ! 258: end ! 259: ! 260: local ! 261: val small=1.0E~9 ! 262: val big=1.0E18 ! 263: val athfhi = 4.6364760900080609352E~1 ! 264: val athflo = 4.6249969567426939759E~18 ! 265: val at1fhi = 9.8279372324732905408E~1 ! 266: val at1flo = ~2.4407677060164810007E~17 ! 267: val a1 = 3.3333333333333942106E~1 ! 268: val a2 = ~1.9999999999979536924E~1 ! 269: val a3 = 1.4285714278004377209E~1 ! 270: val a4 = ~1.1111110579344973814E~1 ! 271: val a5 = 9.0908906105474668324E~2 ! 272: val a6 = ~7.6919217767468239799E~2 ! 273: val a7 = 6.6614695906082474486E~2 ! 274: val a8 = ~5.8358371008508623523E~2 ! 275: val a9 = 4.9850617156082015213E~2 ! 276: val a10 = ~3.6700606902093604877E~2 ! 277: val a11 = 1.6438029044759730479E~2 ! 278: in ! 279: fun atan2(y:real,x:real) = ! 280: let (* copy down the sign of y and x *) ! 281: val signy = copysign(one,y) ! 282: val signx = copysign(one,x) ! 283: exception Done of real ! 284: fun findatan2(x,y,t) = ! 285: let (* begin argument reduction *) ! 286: val k = floor(four * (t+0.0625)) ! 287: val (hi,lo,t) = ! 288: if t < 2.4375 then ! 289: (* truncate 4(t+1/16) to integer for branching *) ! 290: case k of ! 291: (* t is in [0,7/16] *) ! 292: 0 => if t < small ! 293: then raise Done ! 294: (copysign(if signx>zero then t ! 295: else PI-t, ! 296: signy)) ! 297: else (zero,zero,t) ! 298: | 1 => if t < small ! 299: then raise Done ! 300: (copysign(if signx>zero then t ! 301: else PI-t, ! 302: signy)) ! 303: else (zero,zero,t) ! 304: (* t is in [7/16,11/16] *) ! 305: | 2 => (athfhi,athflo,((y+y)-x)/(x+x+y)) ! 306: (* t is in [11/16,19/16] *) ! 307: | 3 => (PIo4,zero,(y-x)/(x+y)) ! 308: | 4 => (PIo4,zero,(y-x)/(x+y)) ! 309: (* t is in [19/16,39/16] *) ! 310: | _ => let val z = y-x ! 311: val y = y+y+y ! 312: val t = x+x ! 313: in (at1fhi,at1flo,( (z+z)-x ) / ( t + y )) ! 314: end ! 315: (* end of if (t < 2.4375) *) ! 316: else let val t = if t <= big ! 317: (* t is in [2.4375, big] *) ! 318: then ~x / y ! 319: (* t is in [big, INF] *) ! 320: else zero ! 321: in (PIo2,zero,t) end ! 322: (* end of argument reduction *) ! 323: (* compute atan(t) for t in [~.4375, .4375] *) ! 324: val z = t*t ! 325: val z = t*(z*(a1+z*(a2+z*(a3+z*(a4+z*(a5+z*(a6+z*(a7+z*(a8+ ! 326: z*(a9+z*(a10+z*a11))))))))))) ! 327: val z = lo - z ! 328: val z = z + t ! 329: val z = z + hi ! 330: in copysign(if signx>zero then z else PI-z,signy) ! 331: end ! 332: in ! 333: (* if x is 1.0, compute *) ! 334: if feql(x,one) then findatan2(x,copysign(y,one),y) handle Done r => r ! 335: (* when y = 0 *) ! 336: else if feql(y,zero) then if feql(signx,one) then y else copysign(PI,signy) ! 337: (* when x = 0 *) ! 338: else if feql(x,zero) then copysign(PIo2,signy) ! 339: else (* compute y/x *) ! 340: let val x = copysign(x,one) ! 341: val y = copysign(y,one) ! 342: val k = logb(y) ! 343: val m = InLine.-(k,logb(x)) ! 344: val (t,y,x) = if InLine.>(m,60) then (big+big,y,x) ! 345: else if InLine.<(m,~80) then (y/x,y,x) ! 346: else (y/x,scalb(y,InLine.~ k), ! 347: scalb(x,InLine.~ k)) ! 348: in findatan2(x,y,t) handle Done r => r end ! 349: end ! 350: end ! 351: ! 352: fun arctan(x:real) = atan2(x,one) ! 353: ! 354: fun sqrt(x: real) = ! 355: if x<=zero then if x<zero then raise Sqrt else x ! 356: else ! 357: let val k = 6 (* log base 2 of the precision *) ! 358: val n = InLine.rshift(logb(x),1) ! 359: val x = scalb(x, InLine.~(InLine.lshift(n,1))) ! 360: fun iter(0,g) = g ! 361: | iter(i,g) = iter(InLine.-(i,1), half * (g + x/g)) ! 362: in scalb(iter(k,one),n) ! 363: end ! 364: ! 365: end
This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.