1*25c28e83SPiotr Jasiukajtis /* 2*25c28e83SPiotr Jasiukajtis * CDDL HEADER START 3*25c28e83SPiotr Jasiukajtis * 4*25c28e83SPiotr Jasiukajtis * The contents of this file are subject to the terms of the 5*25c28e83SPiotr Jasiukajtis * Common Development and Distribution License (the "License"). 6*25c28e83SPiotr Jasiukajtis * You may not use this file except in compliance with the License. 7*25c28e83SPiotr Jasiukajtis * 8*25c28e83SPiotr Jasiukajtis * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE 9*25c28e83SPiotr Jasiukajtis * or http://www.opensolaris.org/os/licensing. 10*25c28e83SPiotr Jasiukajtis * See the License for the specific language governing permissions 11*25c28e83SPiotr Jasiukajtis * and limitations under the License. 12*25c28e83SPiotr Jasiukajtis * 13*25c28e83SPiotr Jasiukajtis * When distributing Covered Code, include this CDDL HEADER in each 14*25c28e83SPiotr Jasiukajtis * file and include the License file at usr/src/OPENSOLARIS.LICENSE. 15*25c28e83SPiotr Jasiukajtis * If applicable, add the following below this CDDL HEADER, with the 16*25c28e83SPiotr Jasiukajtis * fields enclosed by brackets "[]" replaced with your own identifying 17*25c28e83SPiotr Jasiukajtis * information: Portions Copyright [yyyy] [name of copyright owner] 18*25c28e83SPiotr Jasiukajtis * 19*25c28e83SPiotr Jasiukajtis * CDDL HEADER END 20*25c28e83SPiotr Jasiukajtis */ 21*25c28e83SPiotr Jasiukajtis 22*25c28e83SPiotr Jasiukajtis /* 23*25c28e83SPiotr Jasiukajtis * Copyright 2011 Nexenta Systems, Inc. All rights reserved. 24*25c28e83SPiotr Jasiukajtis */ 25*25c28e83SPiotr Jasiukajtis /* 26*25c28e83SPiotr Jasiukajtis * Copyright 2006 Sun Microsystems, Inc. All rights reserved. 27*25c28e83SPiotr Jasiukajtis * Use is subject to license terms. 28*25c28e83SPiotr Jasiukajtis */ 29*25c28e83SPiotr Jasiukajtis 30*25c28e83SPiotr Jasiukajtis #include "libm_synonyms.h" 31*25c28e83SPiotr Jasiukajtis #include "libm_inlines.h" 32*25c28e83SPiotr Jasiukajtis 33*25c28e83SPiotr Jasiukajtis #ifdef __RESTRICT 34*25c28e83SPiotr Jasiukajtis #define restrict _Restrict 35*25c28e83SPiotr Jasiukajtis #else 36*25c28e83SPiotr Jasiukajtis #define restrict 37*25c28e83SPiotr Jasiukajtis #endif 38*25c28e83SPiotr Jasiukajtis 39*25c28e83SPiotr Jasiukajtis #define sqrt __sqrt 40*25c28e83SPiotr Jasiukajtis 41*25c28e83SPiotr Jasiukajtis extern double sqrt(double); 42*25c28e83SPiotr Jasiukajtis 43*25c28e83SPiotr Jasiukajtis void 44*25c28e83SPiotr Jasiukajtis __vhypotf(int n, float * restrict x, int stridex, float * restrict y, 45*25c28e83SPiotr Jasiukajtis int stridey, float * restrict z, int stridez) 46*25c28e83SPiotr Jasiukajtis { 47*25c28e83SPiotr Jasiukajtis float x0, x1, x2, y0, y1, y2, z0, z1, z2, *pz0, *pz1, *pz2; 48*25c28e83SPiotr Jasiukajtis unsigned hx0, hx1, hx2, hy0, hy1, hy2; 49*25c28e83SPiotr Jasiukajtis int i, j0, j1, j2; 50*25c28e83SPiotr Jasiukajtis 51*25c28e83SPiotr Jasiukajtis do 52*25c28e83SPiotr Jasiukajtis { 53*25c28e83SPiotr Jasiukajtis LOOP0: 54*25c28e83SPiotr Jasiukajtis hx0 = *(unsigned*)x & ~0x80000000; 55*25c28e83SPiotr Jasiukajtis hy0 = *(unsigned*)y & ~0x80000000; 56*25c28e83SPiotr Jasiukajtis *(unsigned*)&x0 = hx0; 57*25c28e83SPiotr Jasiukajtis *(unsigned*)&y0 = hy0; 58*25c28e83SPiotr Jasiukajtis if (hy0 > hx0) 59*25c28e83SPiotr Jasiukajtis { 60*25c28e83SPiotr Jasiukajtis i = hy0 - hx0; 61*25c28e83SPiotr Jasiukajtis j0 = hy0 & 0x7f800000; 62*25c28e83SPiotr Jasiukajtis if (hx0 == 0) 63*25c28e83SPiotr Jasiukajtis i = 0x7f800000; 64*25c28e83SPiotr Jasiukajtis } 65*25c28e83SPiotr Jasiukajtis else 66*25c28e83SPiotr Jasiukajtis { 67*25c28e83SPiotr Jasiukajtis i = hx0 - hy0; 68*25c28e83SPiotr Jasiukajtis j0 = hx0 & 0x7f800000; 69*25c28e83SPiotr Jasiukajtis if (hy0 == 0) 70*25c28e83SPiotr Jasiukajtis i = 0x7f800000; 71*25c28e83SPiotr Jasiukajtis else if (hx0 == 0) 72*25c28e83SPiotr Jasiukajtis i = 0x7f800000; 73*25c28e83SPiotr Jasiukajtis } 74*25c28e83SPiotr Jasiukajtis if (i >= 0x0c800000 || j0 >= 0x7f800000) 75*25c28e83SPiotr Jasiukajtis { 76*25c28e83SPiotr Jasiukajtis z0 = x0 + y0; 77*25c28e83SPiotr Jasiukajtis if (hx0 == 0x7f800000) 78*25c28e83SPiotr Jasiukajtis z0 = x0; 79*25c28e83SPiotr Jasiukajtis else if (hy0 == 0x7f800000) 80*25c28e83SPiotr Jasiukajtis z0 = y0; 81*25c28e83SPiotr Jasiukajtis else if (hx0 > 0x7f800000 || hy0 > 0x7f800000) 82*25c28e83SPiotr Jasiukajtis z0 = *x + *y; 83*25c28e83SPiotr Jasiukajtis *z = z0; 84*25c28e83SPiotr Jasiukajtis x += stridex; 85*25c28e83SPiotr Jasiukajtis y += stridey; 86*25c28e83SPiotr Jasiukajtis z += stridez; 87*25c28e83SPiotr Jasiukajtis i = 0; 88*25c28e83SPiotr Jasiukajtis if (--n <= 0) 89*25c28e83SPiotr Jasiukajtis break; 90*25c28e83SPiotr Jasiukajtis goto LOOP0; 91*25c28e83SPiotr Jasiukajtis } 92*25c28e83SPiotr Jasiukajtis pz0 = z; 93*25c28e83SPiotr Jasiukajtis x += stridex; 94*25c28e83SPiotr Jasiukajtis y += stridey; 95*25c28e83SPiotr Jasiukajtis z += stridez; 96*25c28e83SPiotr Jasiukajtis i = 1; 97*25c28e83SPiotr Jasiukajtis if (--n <= 0) 98*25c28e83SPiotr Jasiukajtis break; 99*25c28e83SPiotr Jasiukajtis 100*25c28e83SPiotr Jasiukajtis LOOP1: 101*25c28e83SPiotr Jasiukajtis hx1 = *(unsigned*)x & ~0x80000000; 102*25c28e83SPiotr Jasiukajtis hy1 = *(unsigned*)y & ~0x80000000; 103*25c28e83SPiotr Jasiukajtis *(unsigned*)&x1 = hx1; 104*25c28e83SPiotr Jasiukajtis *(unsigned*)&y1 = hy1; 105*25c28e83SPiotr Jasiukajtis if (hy1 > hx1) 106*25c28e83SPiotr Jasiukajtis { 107*25c28e83SPiotr Jasiukajtis i = hy1 - hx1; 108*25c28e83SPiotr Jasiukajtis j1 = hy1 & 0x7f800000; 109*25c28e83SPiotr Jasiukajtis if (hx1 == 0) 110*25c28e83SPiotr Jasiukajtis i = 0x7f800000; 111*25c28e83SPiotr Jasiukajtis } 112*25c28e83SPiotr Jasiukajtis else 113*25c28e83SPiotr Jasiukajtis { 114*25c28e83SPiotr Jasiukajtis i = hx1 - hy1; 115*25c28e83SPiotr Jasiukajtis j1 = hx1 & 0x7f800000; 116*25c28e83SPiotr Jasiukajtis if (hy1 == 0) 117*25c28e83SPiotr Jasiukajtis i = 0x7f800000; 118*25c28e83SPiotr Jasiukajtis else if (hx1 == 0) 119*25c28e83SPiotr Jasiukajtis i = 0x7f800000; 120*25c28e83SPiotr Jasiukajtis } 121*25c28e83SPiotr Jasiukajtis if (i >= 0x0c800000 || j1 >= 0x7f800000) 122*25c28e83SPiotr Jasiukajtis { 123*25c28e83SPiotr Jasiukajtis z1 = x1 + y1; 124*25c28e83SPiotr Jasiukajtis if (hx1 == 0x7f800000) 125*25c28e83SPiotr Jasiukajtis z1 = x1; 126*25c28e83SPiotr Jasiukajtis else if (hy1 == 0x7f800000) 127*25c28e83SPiotr Jasiukajtis z1 = y1; 128*25c28e83SPiotr Jasiukajtis else if (hx1 > 0x7f800000 || hy1 > 0x7f800000) 129*25c28e83SPiotr Jasiukajtis z1 = *x + *y; 130*25c28e83SPiotr Jasiukajtis *z = z1; 131*25c28e83SPiotr Jasiukajtis x += stridex; 132*25c28e83SPiotr Jasiukajtis y += stridey; 133*25c28e83SPiotr Jasiukajtis z += stridez; 134*25c28e83SPiotr Jasiukajtis i = 1; 135*25c28e83SPiotr Jasiukajtis if (--n <= 0) 136*25c28e83SPiotr Jasiukajtis break; 137*25c28e83SPiotr Jasiukajtis goto LOOP1; 138*25c28e83SPiotr Jasiukajtis } 139*25c28e83SPiotr Jasiukajtis pz1 = z; 140*25c28e83SPiotr Jasiukajtis x += stridex; 141*25c28e83SPiotr Jasiukajtis y += stridey; 142*25c28e83SPiotr Jasiukajtis z += stridez; 143*25c28e83SPiotr Jasiukajtis i = 2; 144*25c28e83SPiotr Jasiukajtis if (--n <= 0) 145*25c28e83SPiotr Jasiukajtis break; 146*25c28e83SPiotr Jasiukajtis 147*25c28e83SPiotr Jasiukajtis LOOP2: 148*25c28e83SPiotr Jasiukajtis hx2 = *(unsigned*)x & ~0x80000000; 149*25c28e83SPiotr Jasiukajtis hy2 = *(unsigned*)y & ~0x80000000; 150*25c28e83SPiotr Jasiukajtis *(unsigned*)&x2 = hx2; 151*25c28e83SPiotr Jasiukajtis *(unsigned*)&y2 = hy2; 152*25c28e83SPiotr Jasiukajtis if (hy2 > hx2) 153*25c28e83SPiotr Jasiukajtis { 154*25c28e83SPiotr Jasiukajtis i = hy2 - hx2; 155*25c28e83SPiotr Jasiukajtis j2 = hy2 & 0x7f800000; 156*25c28e83SPiotr Jasiukajtis if (hx2 == 0) 157*25c28e83SPiotr Jasiukajtis i = 0x7f800000; 158*25c28e83SPiotr Jasiukajtis } 159*25c28e83SPiotr Jasiukajtis else 160*25c28e83SPiotr Jasiukajtis { 161*25c28e83SPiotr Jasiukajtis i = hx2 - hy2; 162*25c28e83SPiotr Jasiukajtis j2 = hx2 & 0x7f800000; 163*25c28e83SPiotr Jasiukajtis if (hy2 == 0) 164*25c28e83SPiotr Jasiukajtis i = 0x7f800000; 165*25c28e83SPiotr Jasiukajtis else if (hx2 == 0) 166*25c28e83SPiotr Jasiukajtis i = 0x7f800000; 167*25c28e83SPiotr Jasiukajtis } 168*25c28e83SPiotr Jasiukajtis if (i >= 0x0c800000 || j2 >= 0x7f800000) 169*25c28e83SPiotr Jasiukajtis { 170*25c28e83SPiotr Jasiukajtis z2 = x2 + y2; 171*25c28e83SPiotr Jasiukajtis if (hx2 == 0x7f800000) 172*25c28e83SPiotr Jasiukajtis z2 = x2; 173*25c28e83SPiotr Jasiukajtis else if (hy2 == 0x7f800000) 174*25c28e83SPiotr Jasiukajtis z2 = y2; 175*25c28e83SPiotr Jasiukajtis else if (hx2 > 0x7f800000 || hy2 > 0x7f800000) 176*25c28e83SPiotr Jasiukajtis z2 = *x + *y; 177*25c28e83SPiotr Jasiukajtis *z = z2; 178*25c28e83SPiotr Jasiukajtis x += stridex; 179*25c28e83SPiotr Jasiukajtis y += stridey; 180*25c28e83SPiotr Jasiukajtis z += stridez; 181*25c28e83SPiotr Jasiukajtis i = 2; 182*25c28e83SPiotr Jasiukajtis if (--n <= 0) 183*25c28e83SPiotr Jasiukajtis break; 184*25c28e83SPiotr Jasiukajtis goto LOOP2; 185*25c28e83SPiotr Jasiukajtis } 186*25c28e83SPiotr Jasiukajtis pz2 = z; 187*25c28e83SPiotr Jasiukajtis 188*25c28e83SPiotr Jasiukajtis z0 = sqrt(x0 * (double)x0 + y0 * (double)y0); 189*25c28e83SPiotr Jasiukajtis z1 = sqrt(x1 * (double)x1 + y1 * (double)y1); 190*25c28e83SPiotr Jasiukajtis z2 = sqrt(x2 * (double)x2 + y2 * (double)y2); 191*25c28e83SPiotr Jasiukajtis *pz0 = z0; 192*25c28e83SPiotr Jasiukajtis *pz1 = z1; 193*25c28e83SPiotr Jasiukajtis *pz2 = z2; 194*25c28e83SPiotr Jasiukajtis 195*25c28e83SPiotr Jasiukajtis x += stridex; 196*25c28e83SPiotr Jasiukajtis y += stridey; 197*25c28e83SPiotr Jasiukajtis z += stridez; 198*25c28e83SPiotr Jasiukajtis i = 0; 199*25c28e83SPiotr Jasiukajtis } while (--n > 0); 200*25c28e83SPiotr Jasiukajtis 201*25c28e83SPiotr Jasiukajtis if (i > 0) 202*25c28e83SPiotr Jasiukajtis { 203*25c28e83SPiotr Jasiukajtis if (i > 1) 204*25c28e83SPiotr Jasiukajtis { 205*25c28e83SPiotr Jasiukajtis z1 = sqrt(x1 * (double)x1 + y1 * (double)y1); 206*25c28e83SPiotr Jasiukajtis *pz1 = z1; 207*25c28e83SPiotr Jasiukajtis } 208*25c28e83SPiotr Jasiukajtis z0 = sqrt(x0 * (double)x0 + y0 * (double)y0); 209*25c28e83SPiotr Jasiukajtis *pz0 = z0; 210*25c28e83SPiotr Jasiukajtis } 211*25c28e83SPiotr Jasiukajtis } 212