PolyMathOrg/PolyMath | Build 4385132063 | src/Math-ArbitraryPrecisionFloat/PMArbitraryPrecisionFloat.class.st | Coveralls

35

PMArbitraryPrecisionFloat class >> example1 [

×

36

×

37

"Compute Pi with nbits of precision"

×

38

×

39

| a b c k pi oldpi oldExpo expo nBits str |

×

40

nBits := 100.

×

41

×

42

a := PMArbitraryPrecisionFloat one asArbitraryPrecisionFloatNumBits: nBits + 16.

×

43

b := (a timesTwoPower: 1) sqrt reciprocal.

×

44

c := a timesTwoPower: -1.

×

45

k := 1.

×

46

oldpi := Float pi.

×

47

oldExpo := 2.

×

48

×

49

[| am gm a2 |

×

50

        am := a + b timesTwoPower: -1.

×

51

        gm := (a * b) sqrt.

×

52

        a := am.

×

53

        b := gm.

×

54

        a2 := a squared.

×

55

        c inPlaceSubtract: (a2 - b squared timesTwoPower: k).

×

56

        pi := (a2 timesTwoPower: 1) / c.

×

57

        expo := (oldpi - pi) exponent.

×

58

        expo isZero or: [expo > oldExpo or: [expo < (-1 - nBits)]]]

×

59

                        whileFalse:

×

60

                                [oldpi := pi.

×

61

                                oldExpo := expo.

×

62

                                k := k + 1].

×

63

pi asArbitraryPrecisionFloatNumBits: nBits.

×

64

str := (String new: 32) writeStream.

×

65

pi absPrintExactlyOn: str base: 10.

×

66

str contents

×

67

×

68

×

69

{ #category : #'instance creation' }

×

70

PMArbitraryPrecisionFloat class >> mantissa: mantisInteger exponent: expoInteger nBits: nbitsInteger [

×

71

        ^self basicNew

×

72

                mantissa: mantisInteger

×

73

                exponent: expoInteger

×

74

                nBits: nbitsInteger

×

75

×

76

×

77

{ #category : #'instance creation' }

×

78

PMArbitraryPrecisionFloat class >> readFrom: aStream numBits: n [

×

79

        "read a number from an ASCII encoded decimal representation with

×

80

        - an optional sign {-|+}

×

81

        - an integer part [0-9]+

×

82

        - an optional decimalPoint and fractionPart {.[0-9]*}

×

83

        - an optional exponent {e{-|+}[0-9]+}"

×

84

×

85

        ^(ExtendedNumberParser on: aStream)

×

86

                failBlock: [self error: 'invalid ArbitraryPrecisionFloat format'];

×

87

                nextArbitraryPrecisionFloatNumBits: n

×

88

×

89

×

90

{ #category : #arithmetic }

×

91

PMArbitraryPrecisionFloat >> * aNumber [

×

92

        | n |

×

93

        aNumber class = self class

×

94

                ifFalse: [^ aNumber adaptToArbitraryPrecisionFloat: self andSend: #*].

×

95

        n := nBits max: aNumber numBits.

×

96

        ^ (self asArbitraryPrecisionFloatNumBits: n)

×

97

                copy inPlaceMultiplyBy: (aNumber asArbitraryPrecisionFloatNumBits: n)

×

98

×

99

×

100

{ #category : #arithmetic }

×

101

PMArbitraryPrecisionFloat >> + aNumber [

×

102

        | n |

×

103

        aNumber class = self class

×

104

                ifFalse: [^ aNumber adaptToArbitraryPrecisionFloat: self andSend: #+].

×

105

        n := nBits max: aNumber numBits.

×

106

        ^ (self asArbitraryPrecisionFloatNumBits: n)

×

107

                copy inPlaceAdd: (aNumber asArbitraryPrecisionFloatNumBits: n)

×

108

×

109

×

110

{ #category : #arithmetic }

×

111

PMArbitraryPrecisionFloat >> - aNumber [

×

112

        | n |

×

113

        aNumber class = self class

×

114

                ifFalse: [^ aNumber adaptToArbitraryPrecisionFloat: self andSend: #-].

×

115

        n := nBits max: aNumber numBits.

×

116

        ^ (self asArbitraryPrecisionFloatNumBits: n)

×

117

                copy inPlaceSubtract: (aNumber asArbitraryPrecisionFloatNumBits: n)

×

118

×

119

×

120

{ #category : #arithmetic }

×

121

PMArbitraryPrecisionFloat >> / aNumber [

×

122

        | n |

×

123

        aNumber class = self class

×

124

                ifFalse: [^ aNumber adaptToArbitraryPrecisionFloat: self andSend: #/].

×

125

        n := nBits max: aNumber numBits.

×

126

        ^ (self asArbitraryPrecisionFloatNumBits: n)

×

127

                copy inPlaceDivideBy: (aNumber asArbitraryPrecisionFloatNumBits: n)

×

128

×

129

×

130

{ #category : #comparing }

×

131

PMArbitraryPrecisionFloat >> < aNumber [

×

132

×

133

        aNumber class = self class ifTrue: [

×

134

                ^ self negative == aNumber negative

×

135

                          ifTrue: [

×

136

                                  self negative

×

137

                                          ifTrue: [ (self digitCompare: aNumber) > 0 ]

×

138

                                          ifFalse: [ (self digitCompare: aNumber) < 0 ] ]

×

139

                          ifFalse: [ self negative ] ].

×

140

        ^ aNumber adaptToArbitraryPrecisionFloat: self andCompare: #<

×

141

×

142

×

143

{ #category : #comparing }

×

144

PMArbitraryPrecisionFloat >> = aNumber [

×

145

×

146

        aNumber isNumber ifFalse: [ ^ false ].

×

147

        aNumber class = self class ifTrue: [

×

148

                ^ aNumber negative == self negative

×

149

                          ifTrue: [ (self digitCompare: aNumber) = 0 ]

×

150

                          ifFalse: [ false ] ].

×

151

        ^ aNumber adaptToArbitraryPrecisionFloat: self andCompare: #=

×

152

×

153

×

154

{ #category : #comparing }

×

155

PMArbitraryPrecisionFloat >> > aNumber [

×

156

×

157

        aNumber class = self class ifTrue: [

×

158

                ^ self negative == aNumber negative

×

159

                          ifTrue: [

×

160

                                  self negative

×

161

                                          ifTrue: [ (self digitCompare: aNumber) < 0 ]

×

162

                                          ifFalse: [ (self digitCompare: aNumber) > 0 ] ]

×

163

                          ifFalse: [ aNumber negative ] ].

×

164

        ^ aNumber adaptToArbitraryPrecisionFloat: self andCompare: #>

×

165

×

166

×

167

{ #category : #printing }

×

168

PMArbitraryPrecisionFloat >> absPrintExactlyOn: aStream base: base [

×

169

        "Print my value on a stream in the given base.

×

170

        Based upon the algorithm outlined in:

×

171

        Robert G. Burger and R. Kent Dybvig

×

172

        Printing Floating Point Numbers Quickly and Accurately

×

173

        ACM SIGPLAN 1996 Conference on Programming Language Design and Implementation

×

174

        June 1996.

×

175

        This version guarantees that the printed representation exactly represents my value

×

176

        by using exact integer arithmetic."

×

177

×

178

        | fBase significand exp baseExpEstimate r s mPlus mMinus scale roundingIncludesLimits d tc1 tc2 fixedFormat decPointCount shead slowbit |

×

179

        fBase := base asFloat.

×

180

        self normalize.

×

181

        significand := mantissa abs.

×

182

        roundingIncludesLimits := significand even.

×

183

        exp := biasedExponent.

×

184

        baseExpEstimate := (self exponent * fBase reciprocalLogBase2 - 1.0e-10) ceiling.

×

185

        exp >= 0

×

186

                ifTrue: [

×

187

                        mPlus := significand isPowerOfTwo

×

188

                                         ifTrue: [

×

189

                                                 r := significand bitShift: 2 + exp.

×

190

                                                 s := 4.

×

191

                                                 2 * (mMinus := 1 bitShift: exp) ]

×

192

                                         ifFalse: [

×

193

                                                 r := significand bitShift: 1 + exp.

×

194

                                                 s := 2.

×

195

                                                 mMinus := 1 bitShift: exp ] ]

×

196

                ifFalse: [

×

197

                        significand isPowerOfTwo

×

198

                                ifTrue: [

×

199

                                        r := significand bitShift: 2.

×

200

                                        s := 1 bitShift: 2 - exp.

×

201

                                        mPlus := 2.

×

202

                                        mMinus := 1 ]

×

203

                                ifFalse: [

×

204

                                        r := significand bitShift: 1.

×

205

                                        s := 1 bitShift: 1 - exp.

×

206

                                        mPlus := mMinus := 1 ] ].

×

207

        baseExpEstimate >= 0

×

208

                ifTrue: [ s := s * (base raisedToInteger: baseExpEstimate) ]

×

209

                ifFalse: [

×

210

                        scale := base raisedToInteger: baseExpEstimate negated.

×

211

                        r := r * scale.

×

212

                        mPlus := mPlus * scale.

×

213

                        mMinus := mMinus * scale ].

×

214

        (r + mPlus >= s and: [ roundingIncludesLimits or: [ r + mPlus > s ] ])

×

215

                ifTrue: [ baseExpEstimate := baseExpEstimate + 1 ]

×

216

                ifFalse: [

×

217

                        r := r * base.

×

218

                        mPlus := mPlus * base.

×

219

                        mMinus := mMinus * base ].

×

220

        (fixedFormat := baseExpEstimate between: -3 and: 6)

×

221

                ifTrue: [

×

222

                        decPointCount := baseExpEstimate.

×

223

                        baseExpEstimate <= 0 ifTrue: [ aStream nextPutAll: ('0.000000' truncateTo: 2 - baseExpEstimate) ] ]

×

224

                ifFalse: [ decPointCount := 1 ].

×

225

        slowbit := 1 - s lowBit.

×

226

        shead := s bitShift: slowbit.

×

227

×

228

        d := (r bitShift: slowbit) // shead.

×

229

        r := r - (d * s).

×

230

        (tc1 := r <= mMinus and: [ roundingIncludesLimits or: [ r < mMinus ] ]) | (tc2 := r + mPlus >= s and: [ roundingIncludesLimits or: [ r + mPlus > s ] ]) ]

×

231

                whileFalse: [

×

232

                        aStream nextPut: (Character digitValue: d).

×

233

                        r := r * base.

×

234

                        mPlus := mPlus * base.

×

235

                        mMinus := mMinus * base.

×

236

                        decPointCount := decPointCount - 1.

×

237

                        decPointCount = 0 ifTrue: [ aStream nextPut: $. ] ].

×

238

        tc2 ifTrue: [ (tc1 not or: [ r * 2 >= s ]) ifTrue: [ d := d + 1 ] ].

×

239

        aStream nextPut: (Character digitValue: d).

×

240

        decPointCount > 0 ifTrue: [

×

241

                decPointCount - 1 to: 1 by: -1 do: [ :i | aStream nextPut: $0 ].

×

242

                aStream nextPutAll: '.0' ].

×

243

        fixedFormat ifFalse: [

×

244

                aStream nextPut: $e.

×

245

                aStream nextPutAll: (baseExpEstimate - 1) printString ]

×

246

×

247

×

248

{ #category : #converting }

×

249

PMArbitraryPrecisionFloat >> adaptToFloat: rcvr andCompare: selector [

×

250

        ^ (rcvr isInfinite or: [ rcvr isNaN ])

×

251

                ifTrue: [ rcvr perform: selector with: self asFloat ]

×

252

                ifFalse: [ (rcvr asArbitraryPrecisionFloatNumBits: 53) perform: selector with: self ]

×

253

×

254

×

255

{ #category : #converting }

×

256

PMArbitraryPrecisionFloat >> adaptToFloat: rcvr andSend: selector [

×

257

        ^ (rcvr isInfinite or: [ rcvr isNaN ])

×

258

                ifTrue: [ rcvr perform: selector with: self asFloat ]

×

259

                ifFalse: [ (rcvr asArbitraryPrecisionFloatNumBits: 53) perform: selector with: self ]

×

260

×

261

×

262

{ #category : #converting }

×

263

PMArbitraryPrecisionFloat >> adaptToFraction: rcvr andCompare: selector [

×

264

        "If I am involved in comparison with a Fraction, convert myself to a

×

265

        Fraction. This way, no bit is lost and comparison is exact."

×

266

×

267

        ^ rcvr perform: selector with: self asTrueFraction

×

268

×

269

×

270

{ #category : #converting }

×

271

PMArbitraryPrecisionFloat >> adaptToFraction: rcvr andSend: selector [

×

272

        ^(rcvr asArbitraryPrecisionFloatNumBits: nBits) perform: selector with: self

×

273

×

274

×

275

{ #category : #converting }

×

276

PMArbitraryPrecisionFloat >> adaptToInteger: rcvr andCompare: selector [

×

277

        "If I am involved in comparison with an Integer, convert myself to a

×

278

        Fraction. This way, no bit is lost and comparison is exact."

×

279

×

280

        ^ rcvr perform: selector with: self asTrueFraction

×

281

×

282

×

283

{ #category : #converting }

×

284

PMArbitraryPrecisionFloat >> adaptToInteger: rcvr andSend: selector [

×

285

        ^(rcvr asArbitraryPrecisionFloatNumBits: nBits) perform: selector with: self

×

286

×

287

×

288

{ #category : #converting }

×

289

PMArbitraryPrecisionFloat >> adaptToScaledDecimal: rcvr andCompare: selector [

×

290

        "If I am involved in comparison with a ScaledDecimal, convert myself to a

×

291

        Fraction. This way, no bit is lost and comparison is exact."

×

292

×

293

        ^ rcvr asFraction perform: selector with: self asFraction

×

294

×

295

×

296

{ #category : #converting }

×

297

PMArbitraryPrecisionFloat >> adaptToScaledDecimal: rcvr andSend: selector [

×

298

        ^(rcvr asArbitraryPrecisionFloatNumBits: nBits) perform: selector with: self

×

299

×

300

×

301

{ #category : #'mathematical functions' }

×

302

PMArbitraryPrecisionFloat >> agm: aNumber [

×

303

        "Answer the arithmetic geometric mean of self and aNumber"

×

304

×

305

        | a b am gm |

×

306

        a := self.

×

307

        b := aNumber.

×

308

×

309

        [am := a + b timesTwoPower: -1.        "am is arithmetic mean"

×

310

        gm := (a * b) sqrt.        "gm is geometric mean"

×

311

        a = am or: [b = gm]]

×

312

                        whileFalse:

×

313

                                [a := am.

×

314

                                b := gm].

×

315

^am

×

316

×

317

×

318

{ #category : #'mathematical functions' }

×

319

PMArbitraryPrecisionFloat >> arCosh [

×

320

        "Evaluate the area hyperbolic cosine of the receiver."

×

321

×

322

        | arCosh x one y two |

×

323

        x := self asArbitraryPrecisionFloatNumBits: 16 + nBits.

×

324

        one := x one.

×

325

        x < one ifTrue: [ DomainError signal: 'cannot compute arCosh of a number less than 1' ].

×

326

        x = one ifTrue: [ ^ self zero ].

×

327

        y := x - one.

×

328

        arCosh := y < one

×

329

                          ifTrue: [

×

330

                                  y exponent * -4 >= nBits

×

331

                                          ifTrue: [ (y powerExpansionArCoshp1Precision: y numBits) * (y timesTwoPower: 1) sqrt ]

×

332

                                          ifFalse: [

×

333

                                                  two := one timesTwoPower: 1.

×

334

                                                  ((y * (y + two)) sqrt + y + one) ln ] ]

×

335

                          ifFalse: [ ((x squared - one) sqrt + x) ln ].

×

336

        ^ arCosh asArbitraryPrecisionFloatNumBits: nBits

×

337

×

338

×

339

{ #category : #'mathematical functions' }

×

340

PMArbitraryPrecisionFloat >> arSinh [

×

341

        "Evaluate the area hyperbolic sine of the receiver."

×

342

×

343

        | arSinh x one |

×

344

        self isZero ifTrue: [ ^ self ].

×

345

        self exponent negated > nBits ifTrue: [ ^ self ].

×

346

        x := self asArbitraryPrecisionFloatNumBits: 16 + nBits.

×

347

        x inPlaceAbs.

×

348

        arSinh := self exponent * -4 >= nBits

×

349

                          ifTrue: [ x powerExpansionArSinhPrecision: x numBits ]

×

350

                          ifFalse: [

×

351

                                  one := x one.

×

352

                                  ((x squared + one) sqrt + x) ln ].

×

353

        self negative ifTrue: [ arSinh inPlaceNegated ].

×

354

        ^ arSinh asArbitraryPrecisionFloatNumBits: nBits

×

355

×

356

×

357

{ #category : #'mathematical functions' }

×

358

PMArbitraryPrecisionFloat >> arTanh [

×

359

        "Evaluate the area hyperbolic tangent of the receiver."

×

360

×

361

        | arTanh x one |

×

362

        self isZero ifTrue: [^self].

×

363

        x := self asArbitraryPrecisionFloatNumBits: 16 + nBits.

×

364

        x inPlaceAbs.

×

365

        one := x one.

×

366

        x >= one ifTrue: [DomainError signal: 'cannot evaluate arTanh of number of magnitude >= 1'].

×

367

        self exponent * -4 >= nBits

×

368

                ifTrue: [arTanh := x powerExpansionArTanhPrecision: x numBits]

×

369

                ifFalse:

×

370

                        [arTanh := ((one + x) / (one - x)) ln.

×

371

                        arTanh inPlaceTimesTwoPower: -1].

×

372

        self negative ifTrue: [arTanh inPlaceNegated].

×

373

        ^arTanh asArbitraryPrecisionFloatNumBits: nBits

×

374

×

375

×

376

{ #category : #'mathematical functions' }

×

377

PMArbitraryPrecisionFloat >> arcCos [

×

378

        "Evaluate the arc cosine of the receiver."

×

379

×

380

        | arcCos x one |

×

381

        self isZero ifTrue: [^(self pi timesTwoPower: -1)].

×

382

        x := self asArbitraryPrecisionFloatNumBits: 16 + nBits.

×

383

        x inPlaceAbs.

×

384

        one := x one.

×

385

        x > one ifTrue: [DomainError signal: 'cannot compute arcCos of a number greater than 1'].

×

386

        arcCos := x = one

×

387

                ifTrue: [self zero]

×

388

                ifFalse: [((one - x squared) sqrt / x) arcTan].

×

389

        self negative ifTrue: [arcCos := x pi - arcCos].

×

390

        ^arcCos asArbitraryPrecisionFloatNumBits: nBits

×

391

×

392

×

393

{ #category : #'mathematical functions' }

×

394

PMArbitraryPrecisionFloat >> arcSin [

×

395

        "Evaluate the arc sine of the receiver."

×

396

×

397

        | arcSin x one |

×

398

        self isZero ifTrue: [^self].

×

399

        x := self asArbitraryPrecisionFloatNumBits: 16 + nBits.

×

400

        x inPlaceAbs.

×

401

        one := x one.

×

402

        x > one ifTrue: [DomainError signal: 'cannot compute arcSin of a number greater than 1'].

×

403

        arcSin := x = one

×

404

                ifTrue: [self pi timesTwoPower: -1]

×

405

                ifFalse: [self exponent * -4 >= nBits

×

406

                        ifTrue: [x powerExpansionArcSinPrecision: x numBits]

×

407

                        ifFalse: [(x / (one - x squared) sqrt) arcTan]].

×

408

        self negative ifTrue: [arcSin inPlaceNegated].

×

409

        ^arcSin asArbitraryPrecisionFloatNumBits: nBits

×

410

×

411

×

412

{ #category : #'mathematical functions' }

×

413

PMArbitraryPrecisionFloat >> arcTan [

×

414

        "Evaluate the arc tangent of the receiver."

×

415

×

416

        | x arcTan one power |

×

417

        self isZero ifTrue: [^self].

×

418

        self > 1

×

419

                ifTrue:

×

420

                        [x := self asArbitraryPrecisionFloatNumBits: nBits * 2 + 2.

×

421

                        x inPlaceAbs.

×

422

                        arcTan := (x pi timesTwoPower: -1) - x reciprocal arcTan]

×

423

                ifFalse:

×

424

                        [power := ((nBits bitShift: -1) + self exponent max: 4) highBit.

×

425

                        x := self asArbitraryPrecisionFloatNumBits: nBits + (1 bitShift: 1 + power).

×

426

                        x inPlaceAbs.

×

427

                        one := x one.

×

428

                        power timesRepeat: [x := x / (one + (one + x squared) sqrt)].

×

429

                        arcTan := x powerExpansionArcTanPrecision: x numBits + 6.

×

430

                        arcTan inPlaceTimesTwoPower: power].

×

431

        self negative ifTrue: [arcTan inPlaceNegated].

×

432

        ^arcTan asArbitraryPrecisionFloatNumBits: nBits

×

433

×

434

×

435

{ #category : #'mathematical functions' }

×

436

PMArbitraryPrecisionFloat >> arcTan: denominator [

×

437

        "Evaluate the four quadrant arc tangent of the argument denominator (x) and the receiver (y)."

×

438

×

439

        ^self isZero

×

440

                ifTrue: [denominator sign positive

×

441

                        ifTrue: [ (self + denominator) zero ]

×

442

                        ifFalse: [ self positive

×

443

                                ifTrue: [ (self + denominator) pi ]

×

444

                                ifFalse: [ (self + denominator) pi negated ]]]

×

445

                ifFalse: [denominator isZero

×

446

                        ifTrue: [self positive

×

447

                                ifTrue: [ (self + denominator) pi timesTwoPower: -1 ]

×

448

                                ifFalse: [ (self + denominator) pi negated timesTwoPower: -1 ]]

×

449

                        ifFalse:

×

450

                                [ | precision arcTan |

×

451

                                precision := (self + denominator) numBits.

×

452

                                arcTan := ((self asArbitraryPrecisionFloatNumBits: precision * 2) / (denominator asArbitraryPrecisionFloatNumBits: precision * 2)) arcTan.

×

453

                                ^(denominator > 0

×

454

                                        ifTrue: [ arcTan ]

×

455

                                        ifFalse: [ self > 0

×

456

                                                ifTrue: [ arcTan + arcTan pi ]

×

457

                                                ifFalse: [ arcTan - arcTan pi ]]) asArbitraryPrecisionFloatNumBits: precision]]

×

458

×

459

×

460

{ #category : #converting }

×

461

PMArbitraryPrecisionFloat >> asArbitraryPrecisionFloatNumBits: n [

×

462

        ^ nBits = n

×

463

                ifTrue: [self]

×

464

                ifFalse: [self copy setPrecisionTo: n]

×

465

×

466

×

467

{ #category : #converting }

×

468

PMArbitraryPrecisionFloat >> asFloat [

×

469

        "Convert to a IEEE 754 double precision floating point."

×

470

×

471

        nBits > Float precision ifTrue: [^(self copy setPrecisionTo: Float precision) asFloat].

×

472

        ^mantissa asFloat timesTwoPower: biasedExponent

×

473

×

474

×

475

{ #category : #converting }

×

476

PMArbitraryPrecisionFloat >> asFraction [

×

477

        ^self asTrueFraction

×

478

×

479

×

480

{ #category : #printing }

×

481

PMArbitraryPrecisionFloat >> asMinimalDecimalFraction [

×

482

        "Answer the shortest decimal Fraction that will equal self when converted back asFloat.

×

483

        A decimal Fraction has only powers of 2 and 5 as decnominator.

×

484

        For example, 0.1 asMinimalDecimalFraction = (1/10)."

×

485

×

486

        | significand exp baseExpEstimate r s mPlus mMinus scale roundingIncludesLimits d tc1 tc2 fixedFormat decPointCount shead slowbit numerator denominator |

×

487

        self isZero ifTrue: [ ^ 0 ].

×

488

        self negative ifTrue: [ ^ self negated asMinimalDecimalFraction negated ].

×

489

        self normalize.

×

490

        significand := mantissa abs.

×

491

        roundingIncludesLimits := significand even.

×

492

        exp := biasedExponent.

×

493

        baseExpEstimate := (self exponent * 10.0 reciprocalLogBase2 - 1.0e-10) ceiling.

×

494

        numerator := 0.

×

495

        denominator := 0.

×

496

        exp >= 0

×

497

                ifTrue: [

×

498

                        mPlus := significand isPowerOfTwo

×

499

                                         ifTrue: [

×

500

                                                 r := significand bitShift: 2 + exp.

×

501

                                                 s := 4.

×

502

                                                 2 * (mMinus := 1 bitShift: exp) ]

×

503

                                         ifFalse: [

×

504

                                                 r := significand bitShift: 1 + exp.

×

505

                                                 s := 2.

×

506

                                                 mMinus := 1 bitShift: exp ] ]

×

507

                ifFalse: [

×

508

                        significand isPowerOfTwo

×

509

                                ifTrue: [

×

510

                                        r := significand bitShift: 2.

×

511

                                        s := 1 bitShift: 2 - exp.

×

512

                                        mPlus := 2.

×

513

                                        mMinus := 1 ]

×

514

                                ifFalse: [

×

515

                                        r := significand bitShift: 1.

×

516

                                        s := 1 bitShift: 1 - exp.

×

517

                                        mPlus := mMinus := 1 ] ].

×

518

        baseExpEstimate >= 0

×

519

                ifTrue: [ s := s * (10 raisedToInteger: baseExpEstimate) ]

×

520

                ifFalse: [

×

521

                        scale := 10 raisedToInteger: baseExpEstimate negated.

×

522

                        r := r * scale.

×

523

                        mPlus := mPlus * scale.

×

524

                        mMinus := mMinus * scale ].

×

525

        (r + mPlus >= s and: [ roundingIncludesLimits or: [ r + mPlus > s ] ])

×

526

                ifTrue: [ baseExpEstimate := baseExpEstimate + 1 ]

×

527

                ifFalse: [

×

528

                        r := r * 10.

×

529

                        mPlus := mPlus * 10.

×

530

                        mMinus := mMinus * 10 ].

×

531

        (fixedFormat := baseExpEstimate between: -3 and: 6)

×

532

                ifTrue: [

×

533

                        decPointCount := baseExpEstimate.

×

534

                        baseExpEstimate <= 0 ifTrue: [ denominator := 10 raisedTo: baseExpEstimate negated ] ]

×

535

                ifFalse: [ decPointCount := 1 ].

×

536

        slowbit := 1 - s lowBit.

×

537

        shead := s bitShift: slowbit.

×

538

×

539

        d := (r bitShift: slowbit) // shead.

×

540

        r := r - (d * s).

×

541

        (tc1 := r <= mMinus and: [ roundingIncludesLimits or: [ r < mMinus ] ]) | (tc2 := r + mPlus >= s and: [ roundingIncludesLimits or: [ r + mPlus > s ] ]) ]

×

542

                whileFalse: [

×

543

                        numerator := 10 * numerator + d.

×

544

                        denominator := 10 * denominator.

×

545

                        r := r * 10.

×

546

                        mPlus := mPlus * 10.

×

547

                        mMinus := mMinus * 10.

×

548

                        decPointCount := decPointCount - 1.

×

549

                        decPointCount = 0 ifTrue: [ denominator := 1 ] ].

×

550

        tc2 ifTrue: [ (tc1 not or: [ r * 2 >= s ]) ifTrue: [ d := d + 1 ] ].

×

551

        numerator := 10 * numerator + d.

×

552

        denominator := 10 * denominator.

×

553

        decPointCount > 0 ifTrue: [ numerator := (10 raisedTo: decPointCount - 1) * numerator ].

×

554

        fixedFormat ifFalse: [

×

555

                baseExpEstimate - 1 > 0

×

556

                        ifTrue: [ numerator := (10 raisedTo: baseExpEstimate - 1) * numerator ]

×

557

                        ifFalse: [ denominator := (10 raisedTo: 1 - baseExpEstimate) * (denominator max: 1) ] ].

×

558

        denominator < 2 ifTrue: [ ^ numerator ].

×

559

        ^ numerator / denominator

×

560

×

561

×

562

{ #category : #converting }

×

563

PMArbitraryPrecisionFloat >> asTrueFraction [

×

564

×

565

        "First remove lowBits from mantissa.

×

566

        This can save a useless division and prevent gcd: cost"

×

567

        self reduce.

×

568

×

569

        ^ biasedExponent >= 0

×

570

                ifTrue: [self shift: mantissa by: biasedExponent]

×

571

                ifFalse: [

×

572

                        "Now avoid a painfull GCD: algorihm.

×

573

                        mantissa is odd and cannot be reduced by a power of two.

×

574

                                mantissa / (1 bitShift: exponent negated)"

×

575

                        ^ Fraction numerator: mantissa denominator: (1 bitShift: biasedExponent negated)]

×

576

×

577

×

578

{ #category : #accessing }

×

579

PMArbitraryPrecisionFloat >> biasedExponent [

×

580

        ^biasedExponent

×

581

×

582

×

583

{ #category : #'mathematical functions' }

×

584

PMArbitraryPrecisionFloat >> cos [

×

585

        ^(PMArbitraryPrecisionFloatForTrigonometry

×

586

                mantissa: mantissa

×

587

                exponent: biasedExponent

×

588

                nBits: nBits) cos

×

589

×

590

×

591

{ #category : #'mathematical functions' }

×

592

PMArbitraryPrecisionFloat >> cosh [

×

593

        | e x |

×

594

        self isZero ifTrue: [^self one].

×

595

        self exponent negated > nBits ifTrue: [^self one].

×

596

        x := self asArbitraryPrecisionFloatNumBits: nBits + 16.

×

597

        self exponent * -4 >= nBits

×

598

                ifTrue: [^(x powerExpansionCoshPrecision: x numBits) asArbitraryPrecisionFloatNumBits: nBits].

×

599

        e := x exp.

×

600

^e

×

601

                inPlaceAdd: e reciprocal;

×

602

                inPlaceTimesTwoPower: -1;

×

603

                asArbitraryPrecisionFloatNumBits: nBits

×

604

×

605

×

606

{ #category : #accessing }

×

607

PMArbitraryPrecisionFloat >> decimalPrecision [

×

608

        ^ nBits * (2 log: 10)

×

609

×

610

×

611

{ #category : #private }

×

612

PMArbitraryPrecisionFloat >> digitCompare: b [

×

613

        "both are positive or negative.

×

614

        answer +1 if i am of greater magnitude, -1 if i am of smaller magnitude, 0 if equal magnitude"

×

615

×

616

        | compare |

×

617

        self isZero ifTrue: [

×

618

                ^ b isZero

×

619

                          ifTrue: [ 0 ]

×

620

                          ifFalse: [ -1 ] ].

×

621

        b isZero ifTrue: [ ^ 1 ].

×

622

        compare := (self exponent - b exponent) sign.

×

623

        ^ compare = 0

×

624

                  ifTrue: [ (self abs - b abs) sign ]

×

625

                  ifFalse: [ compare ]

×

626

×

627

×

628

{ #category : #'mathematical functions' }

×

629

PMArbitraryPrecisionFloat >> exp [

×

630

        "Answer the exponential of the receiver."

×

631

×

632

        | ln2 x q r ri res n maxIter p one two |

×

633

        one := self one.

×

634

        two := one timesTwoPower: 1.

×

635

        "Use following decomposition:

×

636

                x exp = (2 ln * q + r) exp.

×

637

                x exp = (2**q * r exp)"

×

638

        ln2 := two ln.

×

639

        x := self / ln2.

×

640

        q := x truncated.

×

641

        r := (x - q) * ln2.

×

642

×

643

        "now compute r exp by power series expansion

×

644

        we compute (r/(2**p)) exp ** (2**p) in order to have faster convergence"

×

645

        p := 10 min: nBits // 2.

×

646

        r := r timesTwoPower: p negated.

×

647

        ri := one asArbitraryPrecisionFloatNumBits: nBits + 16.

×

648

        res := ri copy.

×

649

        n := 0.

×

650

        maxIter := 1 + ((nBits + 16) / p) ceiling.

×

651

        [n <= maxIter] whileTrue:

×

652

                        [n := n + 1.

×

653

                        ri inPlaceMultiplyBy: r / n.

×

654

                        res inPlaceAdd: ri].

×

655

        p timesRepeat: [res inPlaceMultiplyBy: res].

×

656

        res inPlaceTimesTwoPower: q.

×

657

×

658

        "now use a Newton iteration to refine the result

×

659

        res = res * (self - res ln + 1)"

×

660

        [| oldres delta |

×

661

        oldres := res.

×

662

        res := res asArbitraryPrecisionFloatNumBits: res numBits + 32.

×

663

        res inPlaceMultiplyBy: self - res ln + 1.

×

664

        delta := (res - oldres) exponent.

×

665

        delta = 0 or: [delta <= (res exponent - nBits - 8)]]

×

666

                        whileFalse.

×

667

×

668

        ^res asArbitraryPrecisionFloatNumBits: nBits

×

669

×

670

×

671

{ #category : #accessing }

×

672

PMArbitraryPrecisionFloat >> exponent [

×

673

        "anwser the floating point like exponent e,

×

674

        of self normalized as

×

675

        1.mmmmmm * (2 raisedTo: e)"

×

676

×

677

        self isZero ifTrue: [^0].

×

678

        ^biasedExponent + self numBitsInMantissa - 1

×

679

×

680

×

681

{ #category : #comparing }

×

682

PMArbitraryPrecisionFloat >> hash [

×

683

        "Hash is reimplemented because = is implemented."

×

684

×

685

        ^ self asTrueFraction hash

×

686

×

687

×

688

{ #category : #private }

×

689

PMArbitraryPrecisionFloat >> inPlaceAbs [

×

690

        mantissa := mantissa abs

×

691

×

692

×

693

{ #category : #private }

×

694

PMArbitraryPrecisionFloat >> inPlaceAdd: b [

×

695

        | delta |

×

696

        b isZero ifTrue: [^self round].

×

697

        self isZero

×

698

                ifTrue:

×

699

                        [mantissa := b mantissa.

×

700

                        biasedExponent := b biasedExponent]

×

701

                ifFalse:

×

702

                        [biasedExponent = b biasedExponent

×

703

                                ifTrue: [mantissa := mantissa + b mantissa]

×

704

                                ifFalse:

×

705

                                        ["check for early truncation. beware, keep 2 bits for rounding"

×

706

×

707

                                        delta := self exponent - b exponent.

×

708

                                        delta - 2 > (nBits max: self numBitsInMantissa)

×

709

                                                ifFalse:

×

710

                                                        [delta negated - 2 > (nBits max: b numBitsInMantissa)

×

711

                                                                ifTrue:

×

712

                                                                        [mantissa := b mantissa.

×

713

                                                                        biasedExponent := b biasedExponent]

×

714

                                                                ifFalse:

×

715

                                                                        [delta := biasedExponent - b biasedExponent.

×

716

                                                                        delta > 0

×

717

                                                                                ifTrue:

×

718

                                                                                        [mantissa := (self shift: mantissa by: delta) + b mantissa.

×

719

                                                                                        biasedExponent := biasedExponent - delta]

×

720

                                                                                ifFalse: [mantissa := mantissa + (self shift: b mantissa by: delta negated)]]]]].

×

721

        self round

×

722

×

723

×

724

{ #category : #private }

×

725

PMArbitraryPrecisionFloat >> inPlaceAddNoRound: b [

×

726

        | delta |

×

727

        b isZero ifTrue: [^self].

×

728

        self isZero

×

729

                ifTrue:

×

730

                        [mantissa := b mantissa.

×

731

                        biasedExponent := b biasedExponent]

×

732

                ifFalse:

×

733

                        [delta := biasedExponent - b biasedExponent.

×

734

                        delta isZero

×

735

                                ifTrue: [mantissa := mantissa + b mantissa]

×

736

                                ifFalse:

×

737

                                        [delta > 0

×

738

                                                ifTrue:

×

739

                                                        [mantissa := (self shift: mantissa by: delta) + b mantissa.

×

740

                                                        biasedExponent := biasedExponent - delta]

×

741

                                                ifFalse: [mantissa := mantissa + (self shift: b mantissa by: delta negated)]]]

×

742

×

743

×

744

{ #category : #private }

×

745

PMArbitraryPrecisionFloat >> inPlaceCopy: b [

×

746

        "copy another arbitrary precision float into self"

×

747

×

748

        mantissa := b mantissa.

×

749

        biasedExponent := b biasedExponent.

×

750

        nBits := b numBits

×

751

×

752

×

753

{ #category : #private }

×

754

PMArbitraryPrecisionFloat >> inPlaceDivideBy: y [

×

755

        "Reference: Accelerating Correctly Rounded

×

756

        Floating-Point Division when the Divisor

×

757

        Is Known in Advance - Nicolas Brisebarre,

×

758

        Jean-Michel Muller, Member, IEEE, and

×

759

        Saurabh Kumar Raina -

×

760

        http://perso.ens-lyon.fr/jean-michel.muller/DivIEEETC-aug04.pdf"

×

761

×

762

        | zh x q |

×

763

        zh := y reciprocal reduce.

×

764

        x := self copy.

×

765

        self inPlaceMultiplyBy: zh.

×

766

        q := self copy.

×

767

        "r := "self inPlaceMultiplyBy: y negated andAccumulate: x.

×

768

        "q' := "self inPlaceMultiplyBy: zh andAccumulate: q.

×

769

×

770

        "ALGO 4

×

771

        | zh r zl |

×

772

        zh := b reciprocal.

×

773

        r := b negated inPlaceMultiplyBy: zh andAccumulate: (1 asArbitraryPrecisionFloatNumBits: nBits).

×

774

        zl := (b asArbitraryPrecisionFloatNumBits: nBits + 1) reciprocal inPlaceMultiplyBy: r.

×

775

        self inPlaceMultiplyBy: zh andAccumulate: (zl inPlaceMultiplyBy: self)"

×

776

×

777

×

778

{ #category : #private }

×

779

PMArbitraryPrecisionFloat >> inPlaceMultiplyBy: b [

×

780

        self inPlaceMultiplyNoRoundBy: b.

×

781

        self round

×

782

×

783

×

784

{ #category : #private }

×

785

PMArbitraryPrecisionFloat >> inPlaceMultiplyBy: b andAccumulate: c [

×

786

        "only do rounding after the two operations.

×

787

        This is the traditional muladd operation in aritmetic units"

×

788

×

789

        self inPlaceMultiplyNoRoundBy: b.

×

790

        self inPlaceAdd: c

×

791

×

792

×

793

{ #category : #private }

×

794

PMArbitraryPrecisionFloat >> inPlaceMultiplyNoRoundBy: b [

×

795

        mantissa := mantissa * b mantissa.

×

796

        biasedExponent := biasedExponent + b biasedExponent

×

797

×

798

×

799

{ #category : #private }

×

800

PMArbitraryPrecisionFloat >> inPlaceNegated [

×

801

        mantissa := mantissa negated

×

802

×

803

×

804

{ #category : #private }

×

805

PMArbitraryPrecisionFloat >> inPlaceReciprocal [

×

806

        | ma h |

×

807

        self isZero ifTrue: [(ZeroDivide dividend: self) signal].

×

808

        ma := mantissa abs.

×

809

        h := ma highBit.

×

810

        mantissa := (1 bitShift: h + nBits) + ma quo: (self shift: mantissa by: 1).

×

811

        biasedExponent := biasedExponent negated - h - nBits + 1.

×

812

        self round

×

813

×

814

        "Implementation notes: if m is a power of 2, reciprocal is trivial.

×

815

        Else, we have 2^h > m >2^(h-1)

×

816

        thus 1 < 2^h/m < 2.

×

817

        thus 2^(n-1) < 2^(h+n-1)/m < 2^n

×

818

        We thus have to evaluate (2^(h+n-1)/m) rounded

×

819

        Tie is away from zero because there are always trailing bits (inexact op)

×

820

        (num/den) rounded is also ((num/den)+(sign/2)) truncated

×

821

        or (num*2)+(sign*den) quo: den*2

×

822

        That's finally what we evaluate"

×

823

×

824

×

825

{ #category : #private }

×

826

PMArbitraryPrecisionFloat >> inPlaceSqrt [

×

827

        "Replace the receiver by its square root."

×

828

×

829

        | guess guessSquared delta shift |

×

830

        self negative

×

831

                ifTrue:

×

832

                        [^ DomainError signal: 'sqrt undefined for number less than zero.'].

×

833

        self isZero ifTrue: [^self].

×

834

×

835

        shift := 2 * nBits - mantissa highBit.

×

836

        biasedExponent := biasedExponent - shift.

×

837

        biasedExponent odd

×

838

                ifTrue:

×

839

                        [shift := shift + 1.

×

840

                        biasedExponent := biasedExponent - 1].

×

841

        mantissa := mantissa bitShift: shift.

×

842

        guess := mantissa bitShift: mantissa highBit + 1 // 2.

×

843

×

844

                guessSquared := guess * guess.

×

845

                delta := guessSquared - mantissa quo: (guess bitShift: 1).

×

846

                delta = 0 ] whileFalse:

×

847

                        [ guess := guess - delta ].

×

848

        guessSquared = mantissa

×

849

                ifFalse:

×

850

                        [(guessSquared - guess - mantissa) negative ifFalse: [guess := guess - 1]].

×

851

        mantissa := guess.

×

852

        biasedExponent := biasedExponent quo: 2.

×

853

        self round

×

854

×

855

×

856

{ #category : #private }

×

857

PMArbitraryPrecisionFloat >> inPlaceSubtract: b [

×

858

        | delta |

×

859

        b isZero ifTrue: [^self round].

×

860

        self isZero

×

861

                ifTrue:

×

862

                        [mantissa := b mantissa negated.

×

863

                        biasedExponent := b biasedExponent]

×

864

                ifFalse:

×

865

                        [biasedExponent = b biasedExponent

×

866

                                ifTrue: [mantissa := mantissa - b mantissa]

×

867

                                ifFalse:

×

868

                                        ["check for early truncation. beware, keep 2 bits for rounding"

×

869

×

870

                                        delta := self exponent - b exponent.

×

871

                                        delta - 2 > (nBits max: self numBitsInMantissa)

×

872

                                                ifFalse:

×

873

                                                        [delta negated - 2 > (nBits max: b numBitsInMantissa)

×

874

                                                                ifTrue:

×

875

                                                                        [mantissa := b mantissa negated.

×

876

                                                                        biasedExponent := b biasedExponent]

×

877

                                                                ifFalse:

×

878

                                                                        [delta := biasedExponent - b biasedExponent.

×

879

                                                                        delta >= 0

×

880

                                                                                ifTrue:

×

881

                                                                                        [mantissa := (self shift: mantissa by: delta) - b mantissa.

×

882

                                                                                        biasedExponent := biasedExponent - delta]

×

883

                                                                                ifFalse: [mantissa := mantissa - (self shift: b mantissa by: delta negated)]]]]].

×

884

        self round

×

885

×

886

×

887

{ #category : #private }

×

888

PMArbitraryPrecisionFloat >> inPlaceSubtractNoRound: b [

×

889

        | delta |

×

890

        b isZero ifTrue: [^self].

×

891

        self isZero

×

892

                ifTrue:

×

893

                        [mantissa := b mantissa negated.

×

894

                        biasedExponent := b biasedExponent]

×

895

                ifFalse:

×

896

                        [delta := biasedExponent - b biasedExponent.

×

897

                        delta isZero

×

898

                                ifTrue: [mantissa := mantissa - b mantissa]

×

899

                                ifFalse:

×

900

                                        [delta >= 0

×

901

                                                ifTrue:

×

902

                                                        [mantissa := (self shift: mantissa by: delta) - b mantissa.

×

903

                                                        biasedExponent := biasedExponent - delta]

×

904

                                                ifFalse: [mantissa := mantissa - (self shift: b mantissa by: delta negated)]]]

×

905

×

906

×

907

{ #category : #private }

×

908

PMArbitraryPrecisionFloat >> inPlaceTimesTwoPower: n [

×

909

        self isZero

×

910

                ifFalse: [biasedExponent := biasedExponent + n]

×

911

×

912

×

913

{ #category : #testing }

×

914

PMArbitraryPrecisionFloat >> isAnExactFloat [

×

915

        ^self exponent <= Float emax

×

916

                and: [Float emin - Float precision < self exponent

×

917

                and: [nBits <= Float precision or: [mantissa isAnExactFloat]]]

×

918

×

919

×

920

{ #category : #testing }

×

921

PMArbitraryPrecisionFloat >> isZero [

×

922

        ^mantissa isZero

×

923

×

924

×

925

{ #category : #'mathematical functions' }

×

926

PMArbitraryPrecisionFloat >> ln [

×

927

        "Answer the neperian logarithm of the receiver."

×

928

×

929

        | x4 one two p res selfHighRes prec e |

×

930

        self <= self zero ifTrue: [DomainError signal: 'ln is only defined for x > 0.0'].

×

931

×

932

        one := self one.

×

933

        self = one ifTrue: [^self zero].

×

934

        two := one timesTwoPower: 1.

×

935

×

936

        "Use Salamin algorithm (approximation is good if x is big enough)

×

937

                x ln = Pi  / (2 * (1 agm: 4/x) ).

×

938

        If x not big enough, compute (x timesTwoPower: p) ln - (2 ln * p)

×

939

        if x is close to 1, better use a power expansion"

×

940

        prec := nBits + 16.

×

941

        e := self exponent.

×

942

        e < 0 ifTrue: [e := -1 - e].

×

943

        e > prec

×

944

                ifTrue: [p := 0]

×

945

                ifFalse:

×

946

                        [p := prec - e.

×

947

                        prec := prec + p highBit].

×

948

        selfHighRes := self asArbitraryPrecisionFloatNumBits: prec.

×

949

        (selfHighRes - one) exponent * -4 >= nBits ifTrue: [^(selfHighRes powerExpansionLnPrecision: prec) asArbitraryPrecisionFloatNumBits: nBits].

×

950

        self < one ifTrue: [selfHighRes inPlaceReciprocal].        "Use ln(1/x) => - ln(x)"

×

951

        x4 := (4 asArbitraryPrecisionFloatNumBits: prec)

×

952

                                inPlaceDivideBy: selfHighRes;

×

953

                                inPlaceTimesTwoPower: p negated.

×

954

        res := selfHighRes pi / (one agm: x4) timesTwoPower: -1.

×

955

        res := selfHighRes = two

×

956

                ifTrue: [res / (p + 1)]

×

957

                ifFalse: [p = 0 ifTrue: [res] ifFalse: [res - ((two asArbitraryPrecisionFloatNumBits: prec) ln * p)]].

×

958

        self < one ifTrue: [res inPlaceNegated].

×

959

        ^res asArbitraryPrecisionFloatNumBits: nBits

×

960

×

961

×

962

{ #category : #accessing }

×

963

PMArbitraryPrecisionFloat >> mantissa [

×

964

        ^mantissa

×

965

×

966

×

967

{ #category : #initialization }

×

968

PMArbitraryPrecisionFloat >> mantissa: m exponent: e nBits: n [

×

969

        mantissa := m.

×

970

        biasedExponent := e.

×

971

        nBits := n.

×

972

        self round

×

973

×

974

×

975

{ #category : #private }

×

976

PMArbitraryPrecisionFloat >> moduloNegPiToPi [

×

977

        "answer a copy of the receiver modulo 2*pi, with doubled precision"

×

978

×

979

        | x pi twoPi quo |

×

980

        x := (self asArbitraryPrecisionFloatNumBits: nBits * 2).

×

981

        self negative ifTrue: [x inPlaceNegated].

×

982

        pi := x pi.

×

983

        twoPi := pi timesTwoPower: 1.

×

984

        x > pi ifTrue:

×

985

                [quo := x + pi quo: twoPi.

×

986

                quo highBitOfMagnitude > nBits ifTrue:

×

987

                        [x := (self abs asArbitraryPrecisionFloatNumBits: nBits + quo highBitOfMagnitude).

×

988

                        pi := x pi.

×

989

                        twoPi := pi timesTwoPower: 1.

×

990

                        quo := x + pi quo: twoPi].

×

991

                x inPlaceSubtract: twoPi * quo.

×

992

                self negative ifTrue: [x inPlaceNegated]].

×

993

        ^x asArbitraryPrecisionFloatNumBits: nBits * 2

×

994

×

995

×

996

{ #category : #arithmetic }

×

997

PMArbitraryPrecisionFloat >> naiveRaisedToInteger: n [

×

998

        "Very naive algorithm: use full precision.

×

999

        Use only for small n"

×

1000

        | m e |

×

1001

        m := mantissa raisedToInteger: n.

×

1002

        e := biasedExponent * n.

×

1003

        ^(m asArbitraryPrecisionFloatNumBits: nBits) timesTwoPower: e

×

1004

×

1005

×

1006

{ #category : #arithmetic }

×

1007

PMArbitraryPrecisionFloat >> negated [

×

1008

        ^self copy inPlaceNegated

×

1009

×

1010

×

1011

{ #category : #testing }

×

1012

PMArbitraryPrecisionFloat >> negative [

×

1013

        ^mantissa negative

×

1014

×

1015

×

1016

{ #category : #accessing }

×

1017

PMArbitraryPrecisionFloat >> nextToward: aNumber [

×

1018

        "answer the nearest floating point number to self with same precision than self,

×

1019

        toward the direction of aNumber argument.

×

1020

        If the nearest one falls on the other side of aNumber, than answer a Number"

×

1021

×

1022

        | next |

×

1023

×

1024

        "if self is greater, decrease self, but never under aNumber"

×

1025

        self > aNumber

×

1026

                ifTrue:

×

1027

                        [next := self predecessor.

×

1028

                        ^next >= aNumber

×

1029

                                ifTrue: [next]

×

1030

                                ifFalse: [aNumber]].

×

1031

×

1032

        "if self is smaller, increase self, but never above aNumber"

×

1033

        self < aNumber

×

1034

                ifTrue: [next := self successor.

×

1035

                        ^next <= aNumber

×

1036

                                ifTrue: [next]

×

1037

                                ifFalse: [aNumber]].

×

1038

×

1039

        "if we are equal, return self"

×

1040

        ^self

×

1041

×

1042

×

1043

{ #category : #private }

×

1044

PMArbitraryPrecisionFloat >> normalize [

×

1045

        "normalize the receiver.

×

1046

        a normalized floating point is either 0, or has mantissa highBit = nBits"

×

1047

×

1048

        | delta |

×

1049

        mantissa isZero

×

1050

                ifTrue: [biasedExponent := 0]

×

1051

                ifFalse:

×

1052

                        [self round.

×

1053

                        delta := self numBitsInMantissa - nBits.

×

1054

                        delta < 0

×

1055

                                ifTrue:

×

1056

                                        [mantissa := self shift: mantissa by: delta negated.

×

1057

                                        biasedExponent := biasedExponent + delta]]

×

1058

×

1059

×

1060

{ #category : #accessing }

×

1061

PMArbitraryPrecisionFloat >> numBits [

×

1062

        ^nBits

×

1063

×

1064

×

1065

{ #category : #accessing }

×

1066

PMArbitraryPrecisionFloat >> numBitsInMantissa [

×

1067

        "this is equal to nBits only if we are normalized.

×

1068

        If we are reduced (low bits being zero are removed), then it will be less.

×

1069

        If we haven't been rounded/truncated then it will be more"

×

1070

×

1071

        ^mantissa highBitOfMagnitude

×

1072

×

1073

×

1074

{ #category : #arithmetic }

×

1075

PMArbitraryPrecisionFloat >> one [

×

1076

        ^self class

×

1077

                mantissa: 1

×

1078

                exponent: 0

×

1079

                nBits: nBits

×

1080

×

1081

×

1082

{ #category : #arithmetic }

×

1083

PMArbitraryPrecisionFloat >> pi [

×

1084

        "answer the value of pi rounded to nBits.

×

1085

        Note: use the Brent-Salamin Arithmetic Geometric Mean algorithm"

×

1086

×

1087

        | a b c k pi oldpi oldExpo expo |

×

1088

        a := self one asArbitraryPrecisionFloatNumBits: nBits + 16.

×

1089

        b := (a timesTwoPower: 1) sqrt reciprocal.

×

1090

        c := a timesTwoPower: -1.

×

1091

        k := 1.

×

1092

        oldpi := Float pi.

×

1093

        oldExpo := 2.

×

1094

×

1095

        [| am gm a2 |

×

1096

        am := a + b timesTwoPower: -1.

×

1097

        gm := (a * b) sqrt.

×

1098

        a := am.

×

1099

        b := gm.

×

1100

        a2 := a squared.

×

1101

        c inPlaceSubtract: (a2 - b squared timesTwoPower: k).

×

1102

        pi := (a2 timesTwoPower: 1) / c.

×

1103

        expo := (oldpi - pi) exponent.

×

1104

        expo isZero or: [expo > oldExpo or: [expo < (-1 - nBits)]]]

×

1105

                        whileFalse:

×

1106

                                [oldpi := pi.

×

1107

                                oldExpo := expo.

×

1108

                                k := k + 1].

×

1109

        ^pi asArbitraryPrecisionFloatNumBits: nBits

×

1110

×

1111

×

1112

{ #category : #testing }

×

1113

PMArbitraryPrecisionFloat >> positive [

×

1114

        ^mantissa positive

×

1115

×

1116

×

1117

{ #category : #private }

×

1118

PMArbitraryPrecisionFloat >> powerExpansionArCoshp1Precision: precBits [

×

1119

        "Evaluate arcosh(x+1)/sqrt(2*x) for the receiver x by power series expansion.

×

1120

        The algorithm is interesting when the receiver is close to zero"

×

1121

×

1122

        | one two count count2 sum term term1 term2 |

×

1123

        one := self one.

×

1124

        two := one timesTwoPower: 1.

×

1125

        count := one copy.

×

1126

        count2 := one copy.

×

1127

        sum := one copy.

×

1128

        term1 := one copy.

×

1129

        term2 := one copy.

×

1130

×

1131

        [term1 inPlaceMultiplyBy: self.

×

1132

        term1 inPlaceNegated.

×

1133

        term2 inPlaceMultiplyBy: count2.

×

1134

        term2 inPlaceMultiplyBy: count2.

×

1135

        term2 inPlaceDivideBy: count.

×

1136

        count inPlaceAdd: one.

×

1137

        count2 inPlaceAdd: two.

×

1138

        term2 inPlaceDivideBy: count2.

×

1139

        term2 inPlaceTimesTwoPower: -2.

×

1140

        term := term1 * term2.

×

1141

        sum inPlaceAdd: term.

×

1142

        term exponent + precBits < sum exponent] whileFalse.

×

1143

        ^sum

×

1144

×

1145

×

1146

{ #category : #private }

×

1147

PMArbitraryPrecisionFloat >> powerExpansionArSinhPrecision: precBits [

×

1148

        "Evaluate the area hypebolic sine of the receiver by power series expansion.

×

1149

        The algorithm is interesting when the receiver is close to zero"

×

1150

×

1151

        | one x2 two count sum term |

×

1152

        one := self one.

×

1153

        two := one timesTwoPower: 1.

×

1154

        count := one copy.

×

1155

        sum := one copy.

×

1156

        term := one copy.

×

1157

        x2 := self squared.

×

1158

×

1159

        [term inPlaceMultiplyBy: x2.

×

1160

        term inPlaceMultiplyBy: count.

×

1161

        term inPlaceDivideBy: count + one.

×

1162

        term inPlaceNegated.

×

1163

        count inPlaceAdd: two.

×

1164

        sum inPlaceAdd: term / count.

×

1165

        term exponent + precBits < sum exponent] whileFalse.

×

1166

        sum inPlaceMultiplyBy: self.

×

1167

        ^sum

×

1168

×

1169

×

1170

{ #category : #private }

×

1171

PMArbitraryPrecisionFloat >> powerExpansionArTanhPrecision: precBits [

×

1172

        "Evaluate the area hyperbolic tangent of the receiver by power series expansion.

×

1173

        arTanh (x) = x (1 + x^2/3 + x^4/5 + ... ) for -1 < x < 1

×

1174

        The algorithm is interesting when the receiver is close to zero"

×

1175

×

1176

        | one x2 two count sum term |

×

1177

        one := self one.

×

1178

        two := one timesTwoPower: 1.

×

1179

        count := one copy.

×

1180

        sum := one copy.

×

1181

        term := one copy.

×

1182

        x2 := self squared.

×

1183

×

1184

        [term inPlaceMultiplyBy: x2.

×

1185

        count inPlaceAdd: two.

×

1186

        sum inPlaceAdd: term / count.

×

1187

        term exponent + precBits < sum exponent] whileFalse.

×

1188

        sum inPlaceMultiplyBy: self.

×

1189

        ^sum

×

1190

×

1191

×

1192

{ #category : #private }

×

1193

PMArbitraryPrecisionFloat >> powerExpansionArcSinPrecision: precBits [

×

1194

        "Evaluate the arc sine of the receiver by power series expansion.

×

1195

        The algorithm is interesting when the receiver is close to zero"

×

1196

×

1197

        | one x2 two count sum term |

×

1198

        one := self one.

×

1199

        two := one timesTwoPower: 1.

×

1200

        count := one copy.

×

1201

        sum := one copy.

×

1202

        term := one copy.

×

1203

        x2 := self squared.

×

1204

×

1205

        [term inPlaceMultiplyBy: x2.

×

1206

        term inPlaceMultiplyBy: count.

×

1207

        term inPlaceDivideBy: count + one.

×

1208

        count inPlaceAdd: two.

×

1209

        sum inPlaceAdd: term / count.

×

1210

        term exponent + precBits < sum exponent] whileFalse.

×

1211

        sum inPlaceMultiplyBy: self.

×

1212

        ^sum

×

1213

×

1214

×

1215

{ #category : #private }

×

1216

PMArbitraryPrecisionFloat >> powerExpansionArcTanPrecision: precBits [

×

1217

        "Evaluate the arc tangent of the receiver by power series expansion.

×

1218

        arcTan (x) = x (1 - x^2/3 + x^4/5 - ... ) for -1 < x < 1

×

1219

        The algorithm is interesting when the receiver is close to zero"

×

1220

×

1221

        | count one sum term two x2 |

×

1222

        one := self one.

×

1223

        two := one timesTwoPower: 1.

×

1224

        count := one copy.

×

1225

        sum := one copy.

×

1226

        term := one copy.

×

1227

        x2 := self squared.

×

1228

×

1229

        [term inPlaceMultiplyBy: x2.

×

1230

        term inPlaceNegated.

×

1231

        count inPlaceAdd: two.

×

1232

        sum inPlaceAdd: term / count.

×

1233

        term exponent + precBits < sum exponent] whileFalse.

×

1234

        sum inPlaceMultiplyBy: self.

×

1235

        ^sum

×

1236

×

1237

×

1238

{ #category : #private }

×

1239

PMArbitraryPrecisionFloat >> powerExpansionCosPrecision: precBits [

×

1240

        "Evaluate the cosine of the receiver by power series expansion.

×

1241

        The algorithm is interesting when the receiver is close to zero"

×

1242

×

1243

        | count one sum term two x2 |

×

1244

        one := self one.

×

1245

        two := one timesTwoPower: 1.

×

1246

        count := one copy.

×

1247

        sum := one copy.

×

1248

        term := one copy.

×

1249

        x2 := self squared.

×

1250

×

1251

        [term inPlaceMultiplyBy: x2.

×

1252

        term inPlaceDivideBy: count * (count + one).

×

1253

        term inPlaceNegated.

×

1254

        count inPlaceAdd: two.

×

1255

        sum inPlaceAdd: term.

×

1256

        term exponent + precBits < sum exponent] whileFalse.

×

1257

        ^sum

×

1258

×

1259

×

1260

{ #category : #private }

×

1261

PMArbitraryPrecisionFloat >> powerExpansionCoshPrecision: precBits [

×

1262

        "Evaluate the hyperbolic cosine of the receiver by power series expansion.

×

1263

        The algorithm is interesting when the receiver is close to zero"

×

1264

×

1265

        | count one sum term two x2 |

×

1266

        one := self one.

×

1267

        two := one timesTwoPower: 1.

×

1268

        count := one copy.

×

1269

        sum := one copy.

×

1270

        term := one copy.

×

1271

        x2 := self squared.

×

1272

×

1273

        [term inPlaceMultiplyBy: x2.

×

1274

        term inPlaceDivideBy: count * (count + one).

×

1275

        count inPlaceAdd: two.

×

1276

        sum inPlaceAdd: term.

×

1277

        term exponent + precBits < sum exponent] whileFalse.

×

1278

        ^sum

×

1279

×

1280

×

1281

{ #category : #private }

×

1282

PMArbitraryPrecisionFloat >> powerExpansionLnPrecision: precBits [

×

1283

        "Evaluate the neperian logarithm of the receiver by power series expansion.

×

1284

        For quadratic convergence, use:

×

1285

        ln ((1+y)/(1-y)) = 2 y (1 + y^2/3 + y^4/5 + ... ) = 2 ar tanh( y )

×

1286

        (1+y)/(1-y) = self => y = (self-1)/(self+1)

×

1287

        This algorithm is interesting when the receiver is close to 1"

×

1288

×

1289

        | one |

×

1290

        one := self one.

×

1291

        ^((self - one)/(self + one) powerExpansionArTanhPrecision: precBits) timesTwoPower: 1

×

1292

×

1293

×

1294

{ #category : #private }

×

1295

PMArbitraryPrecisionFloat >> powerExpansionSinCosPrecision: precBits [

×

1296

        "Evaluate the sine and cosine of the receiver by power series expansion.

×

1297

        The algorithm is interesting when the receiver is close to zero"

×

1298

×

1299

        | count one sin cos term |

×

1300

        one := self one.

×

1301

        count := one copy.

×

1302

        cos := one copy.

×

1303

        sin := self copy.

×

1304

        term := self copy.

×

1305

×

1306

        [count inPlaceAdd: one.

×

1307

        term

×

1308

                inPlaceMultiplyBy: self;

×

1309

                inPlaceDivideBy: count;

×

1310

                inPlaceNegated.

×

1311

        cos inPlaceAdd: term.

×

1312

×

1313

        count inPlaceAdd: one.

×

1314

        term

×

1315

                inPlaceMultiplyBy: self;

×

1316

                inPlaceDivideBy: count.

×

1317

        sin inPlaceAdd: term.

×

1318

×

1319

        term exponent + precBits < sin exponent] whileFalse.

×

1320

        ^Array with: sin with: cos

×

1321

×

1322

×

1323

{ #category : #private }

×

1324

PMArbitraryPrecisionFloat >> powerExpansionSinPrecision: precBits [

×

1325

        "Evaluate the sine of the receiver by power series expansion.

×

1326

        The algorithm is interesting when the receiver is close to zero"

×

1327

×

1328

        | count one sum term two x2 |

×

1329

        one := self one.

×

1330

        two := one timesTwoPower: 1.

×

1331

        count := two copy.

×

1332

        sum := self copy.

×

1333

        term := self copy.

×

1334

        x2 := self squared.

×

1335

×

1336

        [term inPlaceMultiplyBy: x2.

×

1337

        term inPlaceDivideBy: count * (count + one).

×

1338

        term inPlaceNegated.

×

1339

        count inPlaceAdd: two.

×

1340

        sum inPlaceAdd: term.

×

1341

        term exponent + precBits < sum exponent] whileFalse.

×

1342

        ^sum

×

1343

×

1344

×

1345

{ #category : #private }

×

1346

PMArbitraryPrecisionFloat >> powerExpansionSinhPrecision: precBits [

×

1347

        "Evaluate the hyperbolic sine of the receiver by power series expansion.

×

1348

        The algorithm is interesting when the receiver is close to zero"

×

1349

×

1350

        | count one sum term two x2 |

×

1351

        one := self one.

×

1352

        two := one timesTwoPower: 1.

×

1353

        count := two copy.

×

1354

        sum := self copy.

×

1355

        term := self copy.

×

1356

        x2 := self squared.

×

1357

×

1358

        [term inPlaceMultiplyBy: x2.

×

1359

        term inPlaceDivideBy: count * (count + one).

×

1360

        count inPlaceAdd: two.

×

1361

        sum inPlaceAdd: term.

×

1362

        term exponent + precBits < sum exponent] whileFalse.

×

1363

        ^sum

×

1364

×

1365

×

1366

{ #category : #private }

×

1367

PMArbitraryPrecisionFloat >> powerExpansionTanPrecision: precBits [

×

1368

        "Evaluate the tangent of the receiver by power series expansion.

×

1369

        The algorithm is interesting when the receiver is close to zero"

×

1370

×

1371

        | count one sum term pow two x2 seidel |

×

1372

        one := self one.

×

1373

        two := one timesTwoPower: 1.

×

1374

        count := two copy.

×

1375

        sum := one copy.

×

1376

        pow := one copy.

×

1377

        x2 := self squared.

×

1378

        seidel := OrderedCollection new: 256.

×

1379

        seidel add: 1.

×

1380

×

1381

        [pow inPlaceMultiplyBy: x2.

×

1382

        pow inPlaceDivideBy: count * (count + one).

×

1383

        count inPlaceAdd: two.

×

1384

        2 to: seidel size do: [:i | seidel at: i put: (seidel at: i-1) + (seidel at: i)].

×

1385

        seidel addLast: seidel last.

×

1386

        seidel size to: 2 by: -1 do: [:i | seidel at: i - 1 put: (seidel at: i-1) + (seidel at: i)].

×

1387

        seidel addFirst: seidel first.

×

1388

        term := pow * seidel first.

×

1389

        sum inPlaceAdd: term.

×

1390

        term exponent + precBits < sum exponent] whileFalse.

×

1391

        sum inPlaceMultiplyBy: self.

×

1392

        ^sum

×

1393

×

1394

×

1395

{ #category : #private }

×

1396

PMArbitraryPrecisionFloat >> powerExpansionTanhPrecision: precBits [

×

1397

        "Evaluate the hyperbolic tangent of the receiver by power series expansion.

×

1398

        The algorithm is interesting when the receiver is close to zero"

×

1399

×

1400

        | count one sum term pow two x2 seidel |

×

1401

        one := self one.

×

1402

        two := one timesTwoPower: 1.

×

1403

        count := two copy.

×

1404

        sum := one copy.

×

1405

        pow := one copy.

×

1406

        x2 := self squared.

×

1407

        seidel := OrderedCollection new: 256.

×

1408

        seidel add: 1.

×

1409

×

1410

        [pow inPlaceMultiplyBy: x2.

×

1411

        pow inPlaceDivideBy: count * (count + one).

×

1412

        pow inPlaceNegated.

×

1413

        count inPlaceAdd: two.

×

1414

        2 to: seidel size do: [:i | seidel at: i put: (seidel at: i-1) + (seidel at: i)].

×

1415

        seidel addLast: seidel last.

×

1416

        seidel size to: 2 by: -1 do: [:i | seidel at: i - 1 put: (seidel at: i-1) + (seidel at: i)].

×

1417

        seidel addFirst: seidel first.

×

1418

        term := pow * seidel first.

×

1419

        sum inPlaceAdd: term.

×

1420

        term exponent + precBits < sum exponent] whileFalse.

×

1421

        sum inPlaceMultiplyBy: self.

×

1422

        ^sum

×

1423

×

1424

×

1425

{ #category : #accessing }

×

1426

PMArbitraryPrecisionFloat >> precision [

×

1427

        ^nBits

×

1428

×

1429

×

1430

{ #category : #'truncation and round off' }

×

1431

PMArbitraryPrecisionFloat >> predecessor [

×

1432

        mantissa = 0 ifTrue: [^self].

×

1433

        mantissa negative ifTrue: [^self negated successor negated].

×

1434

        ^mantissa isPowerOfTwo

×

1435

                ifTrue: [self - (self ulp timesTwoPower: -1)]

×

1436

                ifFalse: [self - self ulp]

×

1437

×

1438

×

1439

{ #category : #printing }

×

1440

PMArbitraryPrecisionFloat >> printOn: aStream [

×

1441

        ^self printOn: aStream base: 10

×

1442

×

1443

×

1444

{ #category : #printing }

×

1445

PMArbitraryPrecisionFloat >> printOn: aStream base: base [

×

1446

        aStream

×

1447

                nextPut: $(;

×

1448

                nextPutAll: self class name;

×

1449

                space;

×

1450

                nextPutAll: #readFrom:;

×

1451

                space;

×

1452

                nextPut: $'.

×

1453

        self negative ifTrue: [aStream nextPut: $-].

×

1454

        self absPrintExactlyOn: aStream base: base.

×

1455

        aStream

×

1456

                nextPut: $';

×

1457

                space;

×

1458

                nextPutAll: #readStream;

×

1459

                space;

×

1460

                nextPutAll: #numBits:;

×

1461

                space;

×

1462

                print: nBits;

×

1463

                nextPut: $)

×

1464

×

1465

×

1466

{ #category : #arithmetic }

×

1467

PMArbitraryPrecisionFloat >> raisedToInteger: anInteger [

×

1468

        | bitProbe highPrecisionSelf n result |

×

1469

        n := anInteger abs.

×

1470

        (n < 5 or: [n * nBits < 512])

×

1471

                ifTrue: [^ self naiveRaisedToInteger: anInteger].

×

1472

        bitProbe := 1 bitShift: n highBit - 1.

×

1473

        highPrecisionSelf := self asArbitraryPrecisionFloatNumBits: n highBit * 2 + nBits + 2.

×

1474

        result := highPrecisionSelf one.

×

1475

×

1476

        [(n bitAnd: bitProbe) = 0 ifFalse: [result := result * highPrecisionSelf].

×

1477

        (bitProbe := bitProbe bitShift: -1) > 0]

×

1478

                whileTrue: [result := result squared].

×

1479

×

1480

        ^ (anInteger negative

×

1481

                ifTrue: [result reciprocal]

×

1482

                ifFalse: [result])

×

1483

                asArbitraryPrecisionFloatNumBits: nBits

×

1484

×

1485

×

1486

{ #category : #arithmetic }

×

1487

PMArbitraryPrecisionFloat >> reciprocal [

×

1488

        ^self copy inPlaceReciprocal

×

1489

×

1490

×

1491

{ #category : #private }

×

1492

PMArbitraryPrecisionFloat >> reduce [

×

1493

        "remove trailing zero bits from mantissa so that we can do arithmetic on smaller integer

×

1494

        (that will un-normalize self)"

×

1495

×

1496

        | trailing |

×

1497

        trailing := mantissa abs lowBit - 1.

×

1498

        trailing > 0

×

1499

                ifFalse: [ ^ self ].

×

1500

        mantissa := self shift: mantissa by: trailing negated.

×

1501

        biasedExponent := biasedExponent + trailing

×

1502

×

1503

×

1504

{ #category : #private }

×

1505

PMArbitraryPrecisionFloat >> round [

×

1506

        "apply algorithm round to nearest even used by IEEE arithmetic"

×

1507

×

1508

        "inexact := ma lowBit <= excess."

×

1509

×

1510

        | excess ma carry |

×

1511

        mantissa isZero

×

1512

                ifTrue: [

×

1513

                        biasedExponent := 0.

×

1514

                        ^ self ].

×

1515

        ma := mantissa abs.

×

1516

        excess := ma highBit - nBits.

×

1517

        excess > 0

×

1518

                ifFalse: [ ^ self ].

×

1519

        carry := ma bitAt: excess.

×

1520

        mantissa := self shift: mantissa by: excess negated.

×

1521

        biasedExponent := biasedExponent + excess.

×

1522

        (carry = 1 and: [ mantissa odd or: [ ma lowBit < excess ] ])

×

1523

                ifFalse: [ ^ self ].

×

1524

        mantissa := mantissa + mantissa sign.

×

1525

        self truncate

×

1526

×

1527

×

1528

{ #category : #initialization }

×

1529

PMArbitraryPrecisionFloat >> setPrecisionTo: n [

×

1530

        nBits := n.

×

1531

        self round

×

1532

×

1533

×

1534

{ #category : #private }

×

1535

PMArbitraryPrecisionFloat >> shift: m by: d [

×

1536

        "shift mantissa m absolute value by some d bits, then restore sign"

×

1537

×

1538

        ^m negative

×

1539

                ifTrue: [(m negated bitShift: d) negated]

×

1540

                ifFalse: [m bitShift: d]

×

1541

×

1542

×

1543

{ #category : #accessing }

×

1544

PMArbitraryPrecisionFloat >> sign [

×

1545

        ^mantissa sign

×

1546

×

1547

×

1548

{ #category : #accessing }

×

1549

PMArbitraryPrecisionFloat >> significand [

×

1550

        ^self timesTwoPower: self exponent negated

×

1551

×

1552

×

1553

{ #category : #accessing }

×

1554

PMArbitraryPrecisionFloat >> significandAsInteger [

×

1555

        self normalize.

×

1556

        ^mantissa abs

×

1557

×

1558

×

1559

{ #category : #'mathematical functions' }

×

1560

PMArbitraryPrecisionFloat >> sin [

×

1561

        ^(PMArbitraryPrecisionFloatForTrigonometry

×

1562

                mantissa: mantissa

×

1563

                exponent: biasedExponent

×

1564

                nBits: nBits) sin

×

1565

×

1566

×

1567

{ #category : #'mathematical functions' }

×

1568

PMArbitraryPrecisionFloat >> sincos [

×

1569

        ^(PMArbitraryPrecisionFloatForTrigonometry

×

1570

                mantissa: mantissa

×

1571

                exponent: biasedExponent

×

1572

                nBits: nBits) sincos

×

1573

×

1574

×

1575

{ #category : #'mathematical functions' }

×

1576

PMArbitraryPrecisionFloat >> sinh [

×

1577

        | e x |

×

1578

        self isZero ifTrue: [^self].

×

1579

        self exponent negated > nBits ifTrue: [^self].

×

1580

        x := self asArbitraryPrecisionFloatNumBits: nBits + 16.

×

1581

        self exponent * -4 >= nBits

×

1582

                ifTrue: [^(x powerExpansionSinhPrecision: x numBits) asArbitraryPrecisionFloatNumBits: nBits].

×

1583

        e := x exp.

×

1584

^e

×

1585

                inPlaceSubtract: e reciprocal;

×

1586

                inPlaceTimesTwoPower: -1;

×

1587

                asArbitraryPrecisionFloatNumBits: nBits

×

1588

×

1589

×

1590

{ #category : #'mathematical functions' }

×

1591

PMArbitraryPrecisionFloat >> sqrt [

×

1592

        "Answer the square root of the receiver."

×

1593

×

1594

        | decimalPlaces n norm guess previousGuess one stopIteration |

×

1595

        self negative

×

1596

                ifTrue:

×

1597

                        [^ DomainError signal: 'sqrt undefined for number less than zero.'].

×

1598

        self isZero ifTrue: [^self].

×

1599

×

1600

        "use additional bits"

×

1601

        decimalPlaces := nBits + 16.

×

1602

        n := self asArbitraryPrecisionFloatNumBits: decimalPlaces.

×

1603

×

1604

        "constants"

×

1605

        one := n one.

×

1606

×

1607

        "normalize n"

×

1608

        norm := n exponent quo: 2.

×

1609

        n := n timesTwoPower: norm * -2.

×

1610

×

1611

        "Initial guess for sqrt(1/n)"

×

1612

        previousGuess := self class

×

1613

                                mantissa: 3

×

1614

                                exponent: -2 - (n exponent quo: 2)

×

1615

                                nBits: decimalPlaces.

×

1616

        guess := previousGuess copy.

×

1617

×

1618

        "use iterations x(k+1) := x*( 1 +  (1-x*x*n)/2) to guess sqrt(1/n)"

×

1619

×

1620

        [guess inPlaceMultiplyNoRoundBy: guess.

×

1621

        guess inPlaceMultiplyBy: n.

×

1622

        guess inPlaceNegated.

×

1623

        guess inPlaceAddNoRound: one.

×

1624

×

1625

        "stop when no evolution of numBits + 12 first bits"

×

1626

        stopIteration := guess isZero or: [guess exponent < (decimalPlaces - 4) negated].

×

1627

        guess inPlaceTimesTwoPower: -1.

×

1628

        guess inPlaceAddNoRound: one.

×

1629

        guess inPlaceMultiplyNoRoundBy: previousGuess.

×

1630

        guess negative ifTrue: [guess inPlaceNegated].

×

1631

×

1632

        guess isZero or: [stopIteration]]

×

1633

                        whileFalse:

×

1634

                                [guess round.

×

1635

                                previousGuess inPlaceCopy: guess].

×

1636

×

1637

        "multiply by n and un-normalize"

×

1638

        guess inPlaceMultiplyBy: n.

×

1639

        guess inPlaceTimesTwoPower: norm.

×

1640

        ^guess asArbitraryPrecisionFloatNumBits: nBits

×

1641

×

1642

×

1643

{ #category : #arithmetic }

×

1644

PMArbitraryPrecisionFloat >> squared [

×

1645

        | result |

×

1646

        result := self copy.

×

1647

        result inPlaceMultiplyBy: self.

×

1648

        ^result

×

1649

×

1650

×

1651

{ #category : #printing }

×

1652

PMArbitraryPrecisionFloat >> storeOn: aStream [

×

1653

        aStream nextPut: $(; nextPutAll: self class name.

×

1654

        aStream space; nextPutAll: 'mantissa:'; space; print: mantissa.

×

1655

        aStream space; nextPutAll: 'exponent:'; space; print: biasedExponent.

×

1656

        aStream space; nextPutAll: 'nBits:'; space; print: nBits.

×

1657

        aStream nextPut: $)

×

1658

×

1659

×

1660

{ #category : #'truncation and round off' }

×

1661

PMArbitraryPrecisionFloat >> successor [

×

1662

        mantissa = 0 ifTrue: [^self].

×

1663

        mantissa negative ifTrue: [^self negated predecessor negated].

×

1664

        ^self + self ulp

×

1665

×

1666

×

1667

{ #category : #'mathematical functions' }

×

1668

PMArbitraryPrecisionFloat >> tan [

×

1669

        ^(PMArbitraryPrecisionFloatForTrigonometry

×

1670

                mantissa: mantissa

×

1671

                exponent: biasedExponent

×

1672

                nBits: nBits) tan

×

1673

×

1674

×

1675

{ #category : #'mathematical functions' }

×

1676

PMArbitraryPrecisionFloat >> tanh [

×

1677

        | e x ep one |

×

1678

        self isZero ifTrue: [^self].

×

1679

        self exponent negated > nBits ifTrue: [^self].

×

1680

        x := self asArbitraryPrecisionFloatNumBits: nBits + 16.

×

1681

        self exponent * -4 >= nBits

×

1682

                ifTrue: [^(x powerExpansionTanhPrecision: x numBits) asArbitraryPrecisionFloatNumBits: nBits].

×

1683

        e := x exp.

×

1684

        one :=x one.

×

1685

        e inPlaceMultiplyBy: e.

×

1686

        ep := e + one.

×

1687

^e

×

1688

                inPlaceSubtract: one;

×

1689

                inPlaceDivideBy: ep;

×

1690

                asArbitraryPrecisionFloatNumBits: nBits

×

1691

×

1692

×

1693

{ #category : #arithmetic }

×

1694

PMArbitraryPrecisionFloat >> timesTwoPower: n [

×

1695

        ^ self isZero

×

1696

                ifTrue: [self]

×

1697

                ifFalse: [self copy inPlaceTimesTwoPower: n]

×

1698

×

1699

×

1700

{ #category : #private }

×

1701

PMArbitraryPrecisionFloat >> truncate [

×

1702

        "remove trailing bits if they exceed our allocated number of bits"

×

1703

×

1704

        | excess |

×

1705

        excess := self numBitsInMantissa - nBits.

×

1706

        excess > 0

×

1707

                ifFalse: [ ^ self ].

×

1708

        mantissa := self shift: mantissa by: excess negated.

×

1709

        biasedExponent := biasedExponent + excess

×

1710

×

1711

×

1712

{ #category : #converting }

×

1713

PMArbitraryPrecisionFloat >> truncated [

×

1714

        "answer the integer that is nearest to self in the interval between zero and self"

×

1715

×

1716

        ^biasedExponent negated > self numBitsInMantissa

×

1717

                ifTrue: [0]

×

1718

                ifFalse: [self shift: mantissa by: biasedExponent]

×

1719

×

1720

×

1721

{ #category : #'truncation and round off' }

×

1722

PMArbitraryPrecisionFloat >> ulp [

×

1723

        mantissa = 0 ifTrue: [^self].

×

1724

        ^self one timesTwoPower: self exponent - nBits + 1

×

1725

×

1726

×

1727

{ #category : #arithmetic }

×

1728

PMArbitraryPrecisionFloat >> zero [

×

1729

        ^self class

×

1730

                mantissa: 0

×

1731

                exponent: 0

×

1732

                nBits: nBits

×

1733

×

PolyMathOrg / PolyMath / 4385132063

Source File
Press 'n' to go to next uncovered line, 'b' for previous

PolyMathOrg / PolyMath / 4385132063

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous