bpo-43420: Simple optimizations for Fraction's arithmetics (GH-24779) · python/cpython@690aca7 (original) (raw)

`@@ -380,32 +380,139 @@ def reverse(b, a):

380

381

`return forward, reverse

382

383

Rational arithmetic algorithms: Knuth, TAOCP, Volume 2, 4.5.1.

384

385

Assume input fractions a and b are normalized.

386

387

1) Consider addition/subtraction.

388

389

Let g = gcd(da, db). Then

390

391

na nb nadb ± nbda

392

a ± b == -- ± -- == ------------- ==

393

da db da*db

394

395

na(db//g) ± nb(da//g) t

396

== ----------------------- == -

397

(da*db)//g d

398

399

Now, if g > 1, we're working with smaller integers.

400

401

Note, that t, (da//g) and (db//g) are pairwise coprime.

402

403

404

removed) and da is coprime with na (since input fractions are

405

normalized), hence (da//g) and na are coprime. By symmetry,

406

(db//g) and nb are coprime too. Then,

407

408

gcd(t, da//g) == gcd(na*(db//g), da//g) == 1

409

gcd(t, db//g) == gcd(nb*(da//g), db//g) == 1

410

411

Above allows us optimize reduction of the result to lowest

412

terms. Indeed,

413

414

g2 = gcd(t, d) == gcd(t, (da//g)(db//g)g) == gcd(t, g)

415

416

t//g2 t//g2

417

a ± b == ----------------------- == ----------------

418

(da//g)(db//g)(g//g2) (da//g)*(db//g2)

419

420

is a normalized fraction. This is useful because the unnormalized

421

denominator d could be much larger than g.

422

423

We should special-case g == 1 (and g2 == 1), since 60.8% of

424

randomly-chosen integers are coprime:

425

https://en.wikipedia.org/wiki/Coprime_integers#Probability_of_coprimality

426

Note, that g2 == 1 always for fractions, obtained from floats: here

427

g is a power of 2 and the unnormalized numerator t is an odd integer.

428

429

2) Consider multiplication

430

431

Let g1 = gcd(na, db) and g2 = gcd(nb, da), then

432

433

nanb nanb (na//g1)*(nb//g2)

434

a*b == ----- == ----- == -----------------

435

dadb dbda (db//g1)*(da//g2)

436

437

Note, that after divisions we're multiplying smaller integers.

438

439

Also, the resulting fraction is normalized, because each of

440

two factors in the numerator is coprime to each of the two factors

441

in the denominator.

442

443

Indeed, pick (na//g1). It's coprime with (da//g2), because input

444

fractions are normalized. It's also coprime with (db//g1), because

445

common factors are removed by g1 == gcd(na, db).

446

447

As for addition/subtraction, we should special-case g1 == 1

448

and g2 == 1 for same reason. That happens also for multiplying

449

rationals, obtained from floats.

450

+

383

451

`def _add(a, b):

384

452

`"""a + b"""

385

da, db = a.denominator, b.denominator

386

return Fraction(a.numerator * db + b.numerator * da,

387

da * db)

453

na, da = a.numerator, a.denominator

454

nb, db = b.numerator, b.denominator

455

g = math.gcd(da, db)

456

if g == 1:

457

return Fraction(na * db + da * nb, da * db, _normalize=False)

458

s = da // g

459

t = na * (db // g) + nb * s

460

g2 = math.gcd(t, g)

461

if g2 == 1:

462

return Fraction(t, s * db, _normalize=False)

463

return Fraction(t // g2, s * (db // g2), _normalize=False)

388

464

389

465

`add, radd = _operator_fallbacks(_add, operator.add)

390

466

391

467

`def _sub(a, b):

392

468

`"""a - b"""

393

da, db = a.denominator, b.denominator

394

return Fraction(a.numerator * db - b.numerator * da,

395

da * db)

469

na, da = a.numerator, a.denominator

470

nb, db = b.numerator, b.denominator

471

g = math.gcd(da, db)

472

if g == 1:

473

return Fraction(na * db - da * nb, da * db, _normalize=False)

474

s = da // g

475

t = na * (db // g) - nb * s

476

g2 = math.gcd(t, g)

477

if g2 == 1:

478

return Fraction(t, s * db, _normalize=False)

479

return Fraction(t // g2, s * (db // g2), _normalize=False)

396

480

397

481

`sub, rsub = _operator_fallbacks(_sub, operator.sub)

398

482

399

483

`def _mul(a, b):

400

484

`"""a * b"""

401

return Fraction(a.numerator * b.numerator, a.denominator * b.denominator)

485

na, da = a.numerator, a.denominator

486

nb, db = b.numerator, b.denominator

487

g1 = math.gcd(na, db)

488

if g1 > 1:

489

na //= g1

490

db //= g1

491

g2 = math.gcd(nb, da)

492

if g2 > 1:

493

nb //= g2

494

da //= g2

495

return Fraction(na * nb, db * da, _normalize=False)

402

496

403

497

`mul, rmul = _operator_fallbacks(_mul, operator.mul)

404

498

405

499

`def _div(a, b):

406

500

`"""a / b"""

407

return Fraction(a.numerator * b.denominator,

408

a.denominator * b.numerator)

501

Same as _mul(), with inversed b.

502

na, da = a.numerator, a.denominator

503

nb, db = b.numerator, b.denominator

504

g1 = math.gcd(na, nb)

505

if g1 > 1:

506

na //= g1

507

nb //= g1

508

g2 = math.gcd(db, da)

509

if g2 > 1:

510

da //= g2

511

db //= g2

512

n, d = na * db, nb * da

513

if d < 0:

514

n, d = -n, -d

515

return Fraction(n, d, _normalize=False)

409

516

410

517

`truediv, rtruediv = _operator_fallbacks(_div, operator.truediv)

411

518

bpo-43420: Simple optimizations for Fraction's arithmetics (GH-24779) · python/cpython@690aca7 (original) (raw)

Rational arithmetic algorithms: Knuth, TAOCP, Volume 2, 4.5.1.

Assume input fractions a and b are normalized.

1) Consider addition/subtraction.

Let g = gcd(da, db). Then

na nb nadb ± nbda

a ± b == -- ± -- == ------------- ==

da db da*db

na*(db//g) ± nb*(da//g) t

== ----------------------- == -

(da*db)//g d

Now, if g > 1, we're working with smaller integers.

Note, that t, (da//g) and (db//g) are pairwise coprime.

Indeed, (da//g) and (db//g) share no common factors (they were

removed) and da is coprime with na (since input fractions are

normalized), hence (da//g) and na are coprime. By symmetry,

(db//g) and nb are coprime too. Then,

gcd(t, da//g) == gcd(na*(db//g), da//g) == 1

gcd(t, db//g) == gcd(nb*(da//g), db//g) == 1

Above allows us optimize reduction of the result to lowest

terms. Indeed,

g2 = gcd(t, d) == gcd(t, (da//g)*(db//g)*g) == gcd(t, g)

t//g2 t//g2

a ± b == ----------------------- == ----------------

(da//g)(db//g)(g//g2) (da//g)*(db//g2)

is a normalized fraction. This is useful because the unnormalized

denominator d could be much larger than g.

We should special-case g == 1 (and g2 == 1), since 60.8% of

randomly-chosen integers are coprime:

https://en.wikipedia.org/wiki/Coprime_integers#Probability_of_coprimality

Note, that g2 == 1 always for fractions, obtained from floats: here

g is a power of 2 and the unnormalized numerator t is an odd integer.

2) Consider multiplication

Let g1 = gcd(na, db) and g2 = gcd(nb, da), then

nanb nanb (na//g1)*(nb//g2)

a*b == ----- == ----- == -----------------

dadb dbda (db//g1)*(da//g2)

Note, that after divisions we're multiplying smaller integers.

Also, the resulting fraction is normalized, because each of

two factors in the numerator is coprime to each of the two factors

in the denominator.

Indeed, pick (na//g1). It's coprime with (da//g2), because input

fractions are normalized. It's also coprime with (db//g1), because

common factors are removed by g1 == gcd(na, db).

As for addition/subtraction, we should special-case g1 == 1

and g2 == 1 for same reason. That happens also for multiplying

rationals, obtained from floats.

Same as _mul(), with inversed b.

na(db//g) ± nb(da//g) t

g2 = gcd(t, d) == gcd(t, (da//g)(db//g)g) == gcd(t, g)