19d5a6349SThomas Gleixner // SPDX-License-Identifier: GPL-2.0-only 21da177e4SLinus Torvalds /* IEEE754 floating point arithmetic 31da177e4SLinus Torvalds * single precision square root 41da177e4SLinus Torvalds */ 51da177e4SLinus Torvalds /* 61da177e4SLinus Torvalds * MIPS floating point support 71da177e4SLinus Torvalds * Copyright (C) 1994-2000 Algorithmics Ltd. 81da177e4SLinus Torvalds */ 91da177e4SLinus Torvalds 101da177e4SLinus Torvalds #include "ieee754sp.h" 111da177e4SLinus Torvalds ieee754sp_sqrt(union ieee754sp x)122209bcb1SRalf Baechleunion ieee754sp ieee754sp_sqrt(union ieee754sp x) 131da177e4SLinus Torvalds { 141da177e4SLinus Torvalds int ix, s, q, m, t, i; 151da177e4SLinus Torvalds unsigned int r; 161da177e4SLinus Torvalds COMPXSP; 171da177e4SLinus Torvalds 181da177e4SLinus Torvalds /* take care of Inf and NaN */ 191da177e4SLinus Torvalds 201da177e4SLinus Torvalds EXPLODEXSP; 219e8bad1fSRalf Baechle ieee754_clearcx(); 221da177e4SLinus Torvalds FLUSHXSP; 231da177e4SLinus Torvalds 241da177e4SLinus Torvalds /* x == INF or NAN? */ 251da177e4SLinus Torvalds switch (xc) { 26d5afa7e9SMaciej W. Rozycki case IEEE754_CLASS_SNAN: 27d5afa7e9SMaciej W. Rozycki return ieee754sp_nanxcpt(x); 28d5afa7e9SMaciej W. Rozycki 291da177e4SLinus Torvalds case IEEE754_CLASS_QNAN: 301da177e4SLinus Torvalds /* sqrt(Nan) = Nan */ 31539bfb57SMaciej W. Rozycki return x; 323f7cac41SRalf Baechle 331da177e4SLinus Torvalds case IEEE754_CLASS_ZERO: 341da177e4SLinus Torvalds /* sqrt(0) = 0 */ 351da177e4SLinus Torvalds return x; 363f7cac41SRalf Baechle 371da177e4SLinus Torvalds case IEEE754_CLASS_INF: 381da177e4SLinus Torvalds if (xs) { 391da177e4SLinus Torvalds /* sqrt(-Inf) = Nan */ 409e8bad1fSRalf Baechle ieee754_setcx(IEEE754_INVALID_OPERATION); 41539bfb57SMaciej W. Rozycki return ieee754sp_indef(); 421da177e4SLinus Torvalds } 431da177e4SLinus Torvalds /* sqrt(+Inf) = Inf */ 441da177e4SLinus Torvalds return x; 453f7cac41SRalf Baechle 461da177e4SLinus Torvalds case IEEE754_CLASS_DNORM: 471da177e4SLinus Torvalds case IEEE754_CLASS_NORM: 481da177e4SLinus Torvalds if (xs) { 491da177e4SLinus Torvalds /* sqrt(-x) = Nan */ 509e8bad1fSRalf Baechle ieee754_setcx(IEEE754_INVALID_OPERATION); 51539bfb57SMaciej W. Rozycki return ieee754sp_indef(); 521da177e4SLinus Torvalds } 531da177e4SLinus Torvalds break; 541da177e4SLinus Torvalds } 551da177e4SLinus Torvalds 561da177e4SLinus Torvalds ix = x.bits; 571da177e4SLinus Torvalds 581da177e4SLinus Torvalds /* normalize x */ 591da177e4SLinus Torvalds m = (ix >> 23); 601da177e4SLinus Torvalds if (m == 0) { /* subnormal x */ 611da177e4SLinus Torvalds for (i = 0; (ix & 0x00800000) == 0; i++) 621da177e4SLinus Torvalds ix <<= 1; 631da177e4SLinus Torvalds m -= i - 1; 641da177e4SLinus Torvalds } 651da177e4SLinus Torvalds m -= 127; /* unbias exponent */ 661da177e4SLinus Torvalds ix = (ix & 0x007fffff) | 0x00800000; 671da177e4SLinus Torvalds if (m & 1) /* odd m, double x to make it even */ 681da177e4SLinus Torvalds ix += ix; 691da177e4SLinus Torvalds m >>= 1; /* m = [m/2] */ 701da177e4SLinus Torvalds 711da177e4SLinus Torvalds /* generate sqrt(x) bit by bit */ 721da177e4SLinus Torvalds ix += ix; 7361100500SAleksandar Markovic s = 0; 7461100500SAleksandar Markovic q = 0; /* q = sqrt(x) */ 751da177e4SLinus Torvalds r = 0x01000000; /* r = moving bit from right to left */ 761da177e4SLinus Torvalds 771da177e4SLinus Torvalds while (r != 0) { 781da177e4SLinus Torvalds t = s + r; 791da177e4SLinus Torvalds if (t <= ix) { 801da177e4SLinus Torvalds s = t + r; 811da177e4SLinus Torvalds ix -= t; 821da177e4SLinus Torvalds q += r; 831da177e4SLinus Torvalds } 841da177e4SLinus Torvalds ix += ix; 851da177e4SLinus Torvalds r >>= 1; 861da177e4SLinus Torvalds } 871da177e4SLinus Torvalds 881da177e4SLinus Torvalds if (ix != 0) { 899e8bad1fSRalf Baechle ieee754_setcx(IEEE754_INEXACT); 901da177e4SLinus Torvalds switch (ieee754_csr.rm) { 9156a64733SRalf Baechle case FPU_CSR_RU: 921da177e4SLinus Torvalds q += 2; 931da177e4SLinus Torvalds break; 9456a64733SRalf Baechle case FPU_CSR_RN: 951da177e4SLinus Torvalds q += (q & 1); 961da177e4SLinus Torvalds break; 971da177e4SLinus Torvalds } 981da177e4SLinus Torvalds } 991da177e4SLinus Torvalds ix = (q >> 1) + 0x3f000000; 1001da177e4SLinus Torvalds ix += (m << 23); 1011da177e4SLinus Torvalds x.bits = ix; 1021da177e4SLinus Torvalds return x; 1031da177e4SLinus Torvalds } 104